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SUMMARY 

Bessel and Hankel functions are widely used in many research and application areas. A 
set of FORTRAN subroutines for calculating these functions is presented. The routines calcu- 
late Bessel and Hankel functions of the first and second kinds, as well as their derivatives, for 
wide ranges of integer order and real or complex argument in single or double precision. 
Depending on the order and argument, one of three evaluation methods is used: the power series 
definition, an Airy function expansion, or an asymptotic expansion. Routines to calculate Airy 
functions and their derivatives are also included. 


1. INTRODUCTION 

Since Bessel and Hankel functions are widely used in many research and application areas, it 
is important to have easy-to-use, reliable computer code to calculate them. The code presented 
here is designed to satisfy that need, in the form of FORTRAN routines to calculate J m , Y m 
(the Bessel functions of first and second kind) and their first derivatives, as well as the Hankel 
functions II ^ and H^ and their derivatives, for wide ranges of order and argument. 

The package consists of a set of FORTRAN routines that calculate Bessel and Hankel func- 
tions with real or complex argument in single or double precision. The user writes a calling 
program that specifies the order m, the argument x, and the function/s to be computed. This 
input is sent to a driver that calls the appropriate routine to do the calculation. 

Three different evaluation methods are used, depending on the combination of input var- 
iables: (I) if the order and argument are sufficient small, the power series definition is used 
directly, (2) if both order and argument are sufficiently large, a series form involving Airy 
function is used, or (3) if the order is small but the argument is large, an asymptotic expansion 
is used. 

An outline of this report follows: Section 2 contains a description of the use of the routines, 
including sample programs and a table of output values. Section 3 lists restrictions on input 
variable ranges. Concluding remarks appear in section 4. The analysis of the three evaluation 
methods appears in appendix A, while appendix B contains the FORTRAN code listing of the 
routines. 


* Retired from NASA Lewis Research Center. 


2. USING THE ROUTINES 


The available routines for calculating Bessel functions in single precision are RBESSEL (real 
argument) and CBESSEL (complex argument); for calculating Hankel functions there are 
RHANKEL (real argument) and CHANKEL (complex argument). The double precision ver- 
sions of these routines are RDBESSEL, CDBESSEL, RDIIANKEL and CDIIANKEL, respec- 
tively. Each routine examines the input variables, decides which evaluation method to use (see 
appendix A), calls an appropriate calculation routine to compute the desired value, and returns 
that value to the calling program through output variables labeled by the user. 

In each case, the inputs are: 

m the integer order of the Bessel/Hankel function 

x the real/complex argument of the Bessel/Hankel function 

icode 1 means compute the Bessel/Hankel function of the first kind 

2 means compute the Bessel/Hankel function of the second kind 

3 means compute the first and second functions and their derivatives 

For each routine, there are four output variables. If icode = I, only the first variable (besj, 
ban I, etc) is calculated; if icode — 2, only the second variable (besy, han2, etc) is calculated; if 
icode = 3, all four variables are calculated. 

In the sample programs that follow, both single and double precision versions of the 
routines are called. It is not necessary to declare both single and double precision variables if 
only one version of the routine is called. However, if icode equals 1 or 2, it is still necessary to 
include all four output variables in the subroutine call. For the sample programs below, note 
that different compilers may require slightly different syntax in places; for example, ‘double 
precision’ may need to be changed to ‘real*8’ and ‘double complex’ to ‘complex* 16’. 

All calculations are done in double precision, but values are returned in either single or 
double precision depending on the calling routine. 


2.1. RBESSEL and RDBESSEL 

Here is a sample program illustrating the use of RBESSEL and RDBESSEL: 

program rbestest 

Integer m, icode 

real besj, besy, dbesj,dbesy,x 

double precision desj,desy,ddesj,ddesy,dx 

x = 1.0 

dx = 1 .OdO 

m = 1 

icode = 3 

call r bessel ( x, m, icode, besj, besy ,d besj, d besy) 
call rdbessel(dx,m, icode, desj,desy,ddesj,ddesy) 
write(6,*) ‘besj = ’,besj,desj 
write(6,*) ‘besy — ’ ,besy,desy 
write(6,*) ‘dbesj = ’,dbesj,ddesj 
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write(6,*) ‘dbesy = ’, dbesy, ddesy 
end 

The output is: 

besj = 0.4400506 0.4400505857449335 

besy = -0.7812128 -0.7812128213002896 

dbesj = 0.3251471 0.3251471008130330 

dbesy = 0.8694698 0.8694697855159653 

2.2. CBESSEL and CDBESSEL 

Here is a sample program illustrating the use of CBESSEL and CDBESSEL: 

program cbestest 
integer in, icode 

complex cbesj,cbesy,cdbesj,cdbesy,x 

double complex cdesj,cdesy,cddesj,cddesy,dx 

x = (1.0, 1.0) 

dx = (1.0d0,1.0d0) 

m = 1 

icode = 3 

call cbessel(x,m, icode, cbesj,cbesy,cdbesj,cdbesy) 

call cdbessel(dx,m, icode, cdesj,cdesy,cddesj,cddesy) 

write(6,*) ‘besj = ’,cbesj,cdesj 

write(6,*) ‘besy = ’,cbesy,cdesy 

write(6,*) ‘dbesj = ’,cdbesj,cddesj 

write(6,*) ‘dbesy = ’,cdbesy,cddesy 

end 

The output is: 

besj = (0.6141603,0.3650281) (0.6141603082191699,0.3650280509978571) 

besy = (-0.6576946,0.6298010) (-0.6576945629749859,0.6298010000344907) 

dbesj = (0.4480143,-0.3719638) (0.4480143229256401,-0.3719637789376947) 

dbesy = (0.4594212, 6.6410817E-02) (0.4594212124350948, 6.6410818437236729E-02) 

2.3. RHANKEL and RDHANKEL 

Here is a sample program illustrating the use of RHANKEL and RDHANKEL: 

program rhantest 
integer m, icode 

complex hanl,han2,dhanl ,dhan2,x 
double complex danl,dan2,ddanl,ddan2,dx 
x = 1.0 
dx = l.OdO 
m — 1 
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icode — 3 

call rhankel(x, m, icode, ban i,han2, dhanl, dhan2) 
call rdhankel (dx,m, icode, dan 1 ,dan2,ddan 1 ,ddan2) 
write(6,*) ‘hanl = ’,hanl,danl 
write(6,*) ‘han2 = 5 ,han2,dan2 
write(6,*) ‘dhanl = ’, dhanl, ddanl 
write(6,*) ‘dhan2 = ’,dhan2,ddan2 
end 

The output is: 

hanl = (0.4400506,-0.7812128) 

ban 2 = (0.4400506,0.7812128) 

dhanl = (0.3251471,0.8694698) 

dhan2 = (0.3251471,-0.8694698) 


(0.4400505857449335,-0.7812128213002896) 

(0.4400505857449335,0.7812128213002896) 

(0.3251471008130330,0.8694697855159653) 

(0.3251471008130330,-0.8694697855159653) 


2.4. CHANKEL and CDHANKEL 

Here is a sample program illustrating the use of CHANKEL and CDIIANKEL: 

program chantest 
integer m, icode 

complex chan 1 ,chaii2,cdhan 1 ,cdhan2,x 
double complex cdanl,cdan2,cddanl,cddan2,dx 
x = (1,0, 1.0) 
dx — (I.0d0,1.0d0) 
m = 1 
icode = 3 

call chankel(x,m, icode, chanl,chan2,cdhanl,cdhan2) 

call cdhankel(dx,m, icode, cdanl,cdan2,cddanl,cddan2) 

write(6,*) ‘hanl = ’,chanl,cdanl 

write(6,*) ‘han2 = ’,chan2,cdan2 

write(6,*) ‘dhanl = ’,cdhanl,cddanl 

write(6,*) ‘dhan2 — ’,cdhan2,cddan2 

end 

The output is: 


hanl = (-1.5640676E-02,-0, 2926665) 

han2 = (1.243961,1.022723) 

dhanl = (0.3816035,8 J457448E-02) 

dhan2 = (0.5144252,-0.8313850) 


(- 1. 5640691815320795 E-02, -0.29266651 19771287) 
(1.243961308253660,1.022722613972843) 
(0,3816035044884034, 8. 7457433497400103E-02) 
(0.5144251413628769,-0.8313849913727896) 


Tables I to IV list Bessel function values at various points. Each of the three evaluation 
methods is represented in the tables. These values were generated on a Unix-based workstation 
using the f77 compiler. Different compilers may give slightly different results. 
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3. CONDITIONS ON INPUT VARIABLES 

The following restrictions apply to the inputs: 

(1) The order m must be a nonnegative integer. 

(2) The argument must be a nonnegative real number x, or a complex number z. 

(3) Avoid parameter values over 32 000 in single precision. 

(4) If x < 2(ERTOLM*m!) A (l/m), the Bessel functions are set to 0 or ERTOLP (machine 
infinity, defined below), as appropriate, due to underflow in the power series and Airy func- 
tion evaluations. For complex argument, replace x by |z|. 

(5) As order and argument both grow, the Airy function calculation breaks down due to 
overflow. The exact region in which this occurs is difficult to describe, but if the returned func- 
tion values are 0.000E f 00 or about ERTOLP (see below), then it should be assumed that 
underflow or overflow has occurred. 

The machine tolerance limits ERTOLP and ERTOLM are set at the beginning of the 
routines RBESSEL, CBESSEL, RDBESSEL and CDBESSEL. ERTOLP and ERTOLM are 
the largest and smallest numbers that don’t cause over- or underflow in single precision. Their 
default values are l.OxlO 35 and l.OxlO" 35 , and can be easily changed in the source code to 
accommodate different computers. 


4. CONCLUSION 

The FORTRAN routines described here calculate the Bessel and Hankel functionsand their 
first derivatives for integer order and real or complex argument insingle or double precision. 
One of three different evaluation methods is used, depending on the order and argument: the 
power series, an Airy functionexpansion, or Hankel’s asymptotic expansion. 

The original single precision, real argument Bessel routines were written by Art Saule. 
These were improved and extended to the other cases by Kevin Kreider. Bruce Clark ran 
numerical tests and suggested further improvements. 
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APPENDIX A 


THE EVALUATION METHODS 
A.l. Choosing the Evaluation Method 

The three evaluation methods (power series, Airy functions, and asymptotic expansions) are 
valid in overlapping regions in the (m,x) plane. The driver routines choose which method to use 
according to the scheme depicted in figures 1 and 2. The difference between the regions in the 
two figures is due to the behavior of the different methods for complex argument off the real 
axis. 

The boundaries for these zones were established by comparing output values of the three 
methods with published results (ref. 1) for a wide range of input values. Since there is consider- 
able overlap in the regions in which the three methods are valid, the final placement of the 
boundaries was guided by judicious interpretation of the data, and is somewhat arbitrary. 


A. 2. Series Truncation 

The three methods used here each involve infinite series, which must betruncated so as to 
provide accuracy to about 14 digits. For programming purposes, these series are put in nested 
form by repeatedly factoring outcommon expressions in the terms. This forces a reverse order 
summation with the final nested term as the starting value, so that the number of terms needed 
to obtain a sufficiently accurate approximation to the series must be known a priori. For each 
series, an index (labeled for each series in sections A.3 to A.5) indicating the proper 
number of terms to be used is developed byexamining the original series: terms are included 
until the next additional term does not change the desired accuracy of the sum. 


A. 3. Evaluation Method 1: Power Series 

The power series expressions for the Bessel functions are useful for calculation if the order 
and argument are both relatively small. Figures 1 and 2 indicate the region in which the power 
series are used for real and complex argument. 

As described in section A.2, all series in the code are written in nested form, and infinite 
series are truncated. The expressions used appear below. 


A. 3.1, Bessel Functions of the First Kind 


From reference 1 (eq. (9.1.10)), 
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A.3.2. Bessel Functions of the Second Kind 

From reference 2 (eq. (77), sec 4.8) (a slightly modified version appears in ref. 1 
(eq. (9.1.11)), 
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L 14 = int(8 + 1.4| z|). 
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(This form of F,^ is chosen since it is convenient for programming using the function routines 
RDFACT and CDFACT.) 

F 0 '(z)=0, f;« - (3-23) 

rz 

A.4. Evaluation Method 2: Airy Function Approximation 

Calculating Bessel functions with relatively large order and argument can be done with 
uniformly asymptotic expressions involving Airy functions. Figures J and 2 indicate the region 
in which these expressions are used for real and complex argument. Calculation of the Airy 
functions is discussed in section 10. 

As described in section A. 2, all series in the code are written in nested form, and infinite 
series are truncated. The expressions used appear below. 
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A.4.1. Bessel Functions of the First Kind 


From reference 1 (eq. (9.3.3)), 
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Several comments must be made here. First, from reference 1 (eqs. (9.3.38) and (9.3.39)), 
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Finally, the coefficients a and b (ref. 1, eq. (9.3.40)) are so small that only a Q , a lt a 2 , b Q , b,, 
and b 2 are needed to obtain the desired accuracy. They are: 
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The approximate expression for J m (mz) is then 
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A. 4. 2. Bessel Functions of the Second Kind 


From reference 1 (eq. (9.2.36)) 
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From reference 1 (eq. (9.3.44)) 
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A. 5. Evaluation Method 3: HankePs Asymptotic Expansion 

HankePs asymptotic expansions can be used to approximate Bessel functions if the order is 
small and the argument large. Figures 1 and 2 indicate the region in which these expansions are 
used for real and complex argument. 

As described in section A.2, all series in the code are written in nested form, and infinite 
series are truncated. The expressions used appear below. 


A.5.1. Bessel Functions of the First Kind 
From reference 1 (eq. (9.2.5)) 
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The derivative is given by reference 1 (eq. (9.2.11)) 


J»« = 


N 


2 

7TZ 


(- R m( z ) sin X - S m( z ) cos x) 


(5.9) 


R in( z ) ~ E I ' 1 ) 1 

k=0 


« 1 - 


4m 2 + 16k 2 - 1 


[4m 2 - (4k - 1; 


(8z) 2 


t - 


(8z) 2 


1 - 


(m,2k) 

(2k) 2k 


1 - 


(5.10) 


J 31 


(8z) 2 


• = 4 ( 4 m 2 + (4k) 2 - l)(4m 2 - (4k - 3) 2 )(4m 2 - (4k - 5) 2 ) k = 1, 2, 
4k (4k - 2) ( 4 m 2 + (4k - 4) 2 - l) 


A k = 


> L 31 
(5.11) 


S m( z ) ~ E (• 1 ) J 

k=0 


4ni 2 + 4(2k + l) 2 - l) 


( 4m 2 - (4k + l) 2 


(rn,2k + 1) 
(2z) 2k+1 


(5.12) 


fa 


4m 2 + 3 


8z 


- b; 

1 - . . . 

b l 

1 - Lsi 



(8z) 2 


(® z ) 2 . 




B k - 


4(4m 2 + (4k + 2) 2 - l)(4m 2 -(4k - 3) 2 ) 

(4m 2 

- (4k - l) 2 ) 

4k(4k + 2) 

( 4 m 2 + (4k - 2) 2 

! - 1 



k = 1, 2, . . . , L 31 

(5.13) 


A.5.2. Bessel Functions of the Second Kind 
From reference 1 (eqs. (9.2.6) and (9.2.12)), 
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[p m (z) c °s X + Q m (z)sin x] 
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[E m (z) c °s X ~ S m (z)sin x] 


A. 6. Airy Functions 

Three different evaluation methods are used to calculate the Airy functions Ai(y) and Bi(y) 
and their derivatives, depending on y = re 1 . For r sufficiently small, the power series is used. 


Otherwise, asymptotic expansions are used - an exponential form if Be 


0 , — 

u 

If, 2, 

3 


3 


, or a 


trigonometric form. Figures 3 and 4 indicate the regions in which these different methods are 
used for real and complex argument. 

As described in section A. 2, all series in the code are written in nested form, and infinite 
series are truncated. The expressions used appear below. 


A.6.1. Power Series 

From reference 1 (eqs. (10.4.2) to (10.4.5)), 

Ai(y) = Cjf(y) - c 2 g(y) 


Bi(y) = / 3* (c 1 f(y) + c 2 g(y)) 
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c 2 = 0.25881 94037 92807 


( 6 . 1 ) 

( 6 . 2 ) 

(6.3) 

(6.4) 


f (y) = £ 3 1 

k=Q 


( \ 
1 


3k 


(3k)! 

w 1 + Fjy 3 


1 + F 2 y 3 


1 + 


1 + F L..y 


26 


(6.5) 


F k = 


3k(3k - 1) 


k = 1, 2, . . . , L 


25 


(6.6) 


16 


g(y) = E 3 * 

k=0 


« y 


r 3k+l 


Uik (3k + 1)! 


1 + G x y : 


1 + G 2 y 


3 [l + . . . 


1 + G L„,y 


'26 


(6.7) 


G k = 


3k(3k + 1) 


k = 1, 2, , L 


25 
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A 


( 6 . 10 ) 


The power series has limited usefulness in the complex domain due to cancellation error 
caused by the subtraction in (6.1), hence the limiting radius is lower than for strictly real 
argument. 

For the derivatives, 


Ai'(y) = Cl l'(y) - c 2 g'(y) 


( 6 . 11 ) 


Bi'(y) = (cji'ty) + c 2 g'(y)) 
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A. 6 . 2 . Trigonometric Asymptotic Expansion 
From reference 1 (eqs. (10.4.60) and (10.4.64)) 
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(6.15) 
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For the derivatives, use reference 1 (eqs. (10.4.62) and (10.4,67)), 
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A.6.3. Exponential Approximation 
From reference 1 (eqs. (10.4.59) and (10.4.63)), 
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For the derivatives, use reference 1 (eqs. (10.4.61) and (10.4.66)), 
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Appendix B. FORTRAN Code Listing 

Listed below is the FORTRAN code listing of all the routines. Related routines are grouped 
together; the single precision, real argument routines appear First, followed by the single 
precision, complex argument routines, and then the double precision real and complex versions. 


Single Precision, Real Argument Routines 


RBESSEL 

RHANKEL 

RBSJl 

RDBSJ1 

RBSYl 

RDBSYl 

RBESJY2 

RBESJY3 


Bessel function driver, user-called 
Hankel function driver, user-called 
Bessel function J m (x) via power series 

Bessel function J^,(x) via power series 
Bessel function Y m (x) via power series 

Bessel function Y ni (x) via power series 

all Bessel functions via Airy function approximation 

all Bessel functions via asymptotic expansion 


Single Precision, Complex Argument Routines 


CBESSEL 

CIIANKEL 

CBSJl 


CDBSJ1 

CBSYl 

CDBSY1 

CBESJY2 

CBESJY3 


Bessel function driver, user-called 
Hankel function driver, user-called 
Bessel function J m (x) via power series 

Bessel function J in (x) via power series 
Bessel function Y m (x) via power series 

Bessel function Y m (x) via power series 1 : - 1 

all Bessel functions via Airy function approximation 
all Bessel functions via asymptotic expansion 


Double Precision, Real Argument Routines 


RDBESSEL 
RDHANKEL 
D RBSJl 
RDFACT 

DRDBSJl 

DRBSYl 

DDIGAM 

DRDBSY1 

DRBESJY2 

RDAIRYFN 

DRBESJY3 

RDTOLCH 


Bessel function driver, user-called 
Hankel function driver, user-called 
Bessel function J (x) via power series 
utility: (x/2)‘n/n! 

Bessel function J m (x) via power series 
Bessel function Y ni (x) via power series 
utility: digamma function 

Bessel function Y^x) via power series 

all Bessel functions via Airy function approximation 

utility: Airy functions 

all Bessel functions via asymptotic expansion 
utility: check for under-overflow 
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Double Precision, Complex Argument Routines 


CDBESSEL 
CDHANKEL 
DCBSJ 1 
CDFACT 

DCDBSJ1 

DCBSY1 

DDIGAM 

DCDBSY1 

DCBESJY2 

CDAIRYFN 

DCBESJY3 

CDTOLCII 


Bessel function driver, user-called 
Hankel function driver, user-called 
Bessel function J m (x) via power series 
utility: (x/2)‘n/n! 

Bessel function J^x) via power series 
Bessel function Y m (x) via power series 
utility: digamma function 

Bessel function Y^x) via power series 

all Bessel functions via Airy function approximation 

utility: Airy functions 

all Bessel functions via asymptotic expansion 
utility: check for under-overflow 


*************** *************************** 

subroutine rbessel(arg,m, icode, besj, besy ,dbesj,dbesy) 

c This set of routines calculates various types of Bessel functions 
c in single precision. 

c INPUTS: 

c arg — real argument x >= 0 

c m — integer order m >= 0 

c icode — 1 calculate Jm(x), the Bessel fen of first kind 

c 2 Ym(x), the Bessel fen of second kiud 

c 3 Jm(x), Ym(x) and their First derivatives 

c OUTPUT: 

c the value of the indicated Bessel function, written to the following 
c variables, according to the value of icode: 
c icode = 1 — > besj 

c 2 besy 

c 3 besj, besy ,dbesj,dbesy 

c NOTE: Ym(0) and Y’m(O) are undefined. 

c Three different evaluation methods are used, depending on the 
c relationship between M and ARG: 

c 1) if ARG is small, the power series definition is used directly 
c (routines RBSJl, RDBSJl, RBSY1, RDBSYl). 

c 2) if both are large, a series involving Airy functions is used 
c (routine RBESJY3). 

c 3) if M is small but ARG is large, an asymptotic expansion is used 
c (routine RBESJY2). 
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Integer m, icode 

real arg, zarg, zm,besj,dbesj,besy,dbesy 
double precision ertolp,ertolm 

common /DERRTOL/ ertolp,ertolm 

ertolp = 1.0e35 
ertolm = 1.0e-35 

\ - - 

zarg = min( .007*m**2+.16*m +12.5, 40. ) 
zm = min( 3.*arg/ll.+20., 33. ) 

if (arg .le. zarg) then 
c use power series 

if (icode .eq. l) call rbsjl(arg,m,besj) 
if (icode .eq. 2) call rbsy i(arg,m,besy) 
if (icode .eq. 3) then 
call rbsjl (arg, m,besj) 
call rdbsjl(arg,m,dbesj) 
call rbsy i (arg, m,besy) 
call rdbsyl(arg,m,dbesy) 
endif 

else if (m .ge. zm) then 
c use airy function approximation 

call rbesjy2(arg,m, icode, besj,besy,dbesj,dbesy) 
else 

c use asymptotic expansion approximation 

call rbesjy 3(arg, m, icode, besj,besy,dbesj,dbesy) 
endif 

return 

end 

subroutine rhankel(arg, order, icode, hanl,han2,dhanl,dhan2) 

c This set of routines calculates various types of Hankel functions 
c in single precision. 

c INPUTS: 

c arg — real argument x >= 0 

c order -- integer order m >= 0 

c icode — 1 calculate H[l]rn(x), the Hankel fen of the first kind 

c 2 H[2|m(x), the Hankel fen of the second kind 

c 3 H[l]m(x), H[2]m(x) & first derivatives 

c OUTPUT: 

c the value of the indicated Hankel function, written to the following 
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c variables, according to the value of icode: 
c icode = 1 --> haul 

c 2 han2 

c 3 Iianl,han2,dhanl,dhan2 

c NOTE: the variables hanl,han2,dhanl ,dhan2 must be declared complex 
c in the calling program. 

integer order, icode 

complex han 1 ,dhan 1 ,han2,dhan2,ic 

real arg,besj,dbesj,besy,dbesy 

ic = (0.,1.) 

if (arg .ge. 0.) then 

call rbessel(arg,order,3,besj,besy,dbesj,dbesy) 
if (icode .eq. 1) then 
hanl = besj + ic*besy 
else if (icode .eq. 2) then 
han2 — besj - ic*besy 
else 

hanl = besj + ic*besy 
han2 = besj - ic*besy 
dhanl = dbesj + ic*dbesy 
dhan2 = dbesj - ic*dbesy 
endif 
else 

arg = -arg 

c Abr-Stegun 9.1.39 

call rbessel(arg, order, 3, besj, besy, dbesj, dbesy) 
if (icode .eq. 1) then 

hanl = -(-l.)**order*(besj-ic*besy) 
else if (icode .eq. 2) then 

han2 = -(-l.)**order*(besj+ic*besy) 
else 

hanl = -(-l.)**order*(besj-ic*besy) 
han2 = -(-L)**order*(besj+ic*besy) 
dhanl = (-l.)**order*(dbesj-ic*dbesy) 
dhari2 = (-l.)**order*(dbesj+ic*dbesy) 
endif 
endif 

return 

end 

*** ***************************** * ******** * 
subroutine rbsj 1 (arg, m, besj 1) 

integer m 

double precision darg,besj,ertolp,ertolm 
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real arg,besjl 

common /DERRTOL/ ertolp,ertolm 
darg = arg 

call drbsji(darg,m,besj) 
besjl = besj 

return 

end 

^ 4 : 4 : 4 ! 4 : 4 : 4 :* 4 : 4 :***** + ******************** : 4 : *** + *** 

subroutine rdbsjl(arg,m,dbesjl) 
integer m 

double precision darg,dbesj,ertolp,ertolm 
real arg,dbesjl 

common /DERRTOL/ ertolp,ertolm 
darg — arg 

call drdbsjl(darg,m,dbesj) 
d besjl = dbesj 

return 

end 

subroutine rbsyJ (arg,m,besyl) 
integer m 

double precision darg,besy,ertolp,ertolm 
real arg,besyl 

common /DERRTOL/ erto!p,ertolm 
darg = arg 

call drbsy l(darg,rn,besy) 
besyl = besy 

return 

end 

* + 4: 4< 4 4c 4 ; * 4 : ^ * * * * * 4: 4 * * 4: * * * + 4: 4: + * 4: * + * * * * * * * * * * * 

subroutine rdbsy l(arg,m,dbesy 1) 
integer m 

double precision darg,dbesy,ertolp,ertolm 
real arg,dbesyl 
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common /DERRTOL/ ertolp,ertolin 


darg = arg 

call drdbsyl(darg,m, dbesy) 
dbesyl = dbesy 

return 

end 

subroutine rbesjy2(arg,m,icode,besj2,besy2,dbesj2,dbesy2) 
c Airy function approximation to Bessel function, real argument 
integer m,icode 

double precision darg, besj,besy,dbesj, dbesy ,ertolp,ertolm 
real arg, besj 2 ,besy 2 ,dbesj2 ,d besy 2 

common /DERRTOL/ ertolp,ertolm 

darg — arg 

call drbesjy2(darg,m,icode, besj, besy ,dbesj, dbesy) 

besj 2 = besj 

besy 2 = besy 

dbesj2 = dbesj 

dbesy 2 = dbesy 

return 

end 

subroutine rbesjy3(arg,m,icode,besj3,besy3,dbesj3,dbesy3) 
c asymptotic approximation to Bessel function, real argument 
integer m,icode 

double precision darg, besj, besy, dbesj, dbesy ,ertolp,ertolm 
real arg, besj 3, besy 3, dbesj3, dbesy 3 

common /DERRTOL/ ertolp,ertolm 

darg = arg 

call dr besj y 3 (dar g,m,icode,besj ,besy , dbesj , dbesy ) 

besj 3 = besj 

besy 3 = besy 

dbesj 3 = dbesj 

dbesy3 — dbesy 

return 

end 
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c 

subroutine cbessel(arg, m, icode, besjc, besyc, dbesjc,dbesyc) 

c This set of routines calculates various types of Bessel functions 
c in single precision. 

c INPUTS: 

c arg — complex argument 

cm— integer order m >= 0 

c icode — 1 calculate Jm(x), the Bessel fen of first kind 

c 2 Ym(x), the Bessel fen of second kind 

c 3 Jm(x),Ym(x) & their first derivatives 

c OUTPUT: 

c the value of the indicated Bessel function, written to the following 
c variables, according to the value of icode: 
c icode = 1 — > besjc 

c 2 besyc 

c 3 besjc, besyc, dbesjc,dbesyc 

c NOTE: Ym(0) and Y’m(O) are undefined. 

c Three different evaluation methods are used, depending on the 
c relationship between M and ARG: 

c 1) if both are small, the power series definition is used directly 
c (routines CBSJl, CDBSJ1, CBSY1, CDBSYl). 

c 2) if both are large, a series involving Airy functions is used 
c (routine CBESJY3). 

c 3) if M is small but ARG is large, an asymptotic expansion is used 
c (routine CBESJY2). 

integer m, icode 

complex arg, besjc, besyc, dbesjc,dbesyc 
real zm 

double precision ertolp,ertolm 

common /DERRTOL/ ertolp,ertolm 

ertolp = 1.0d35 
ertolm — 1.0d-35 

if (m .It. 8) then 

zm = 7. + 2.*float(rn)/7. 
else 

zm = (184. - 8.*float(m))/15. 
endif 

if (abs(arg) .le. zm) then 
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c 


c 


c 


use power series 

if (icode .eq. 1) call cbsjl(arg,m,besjc) 
if (icode .eq. 2) call cbsy l(arg,m,besyc) 
if (icode .eq. 3) then 

call cbsjl(arg,m,besjc) 
call cdbsj 1 (arg,m,dbesjc) 
call cbsy 1 (arg, m,besyc) 
call cdbsy l(arg,m,dbesyc) 

endif 

else if (m .ge. 8) then 
use airy approximation 

call cbesjy2(arg,m,icode,besjc,besyc,dbesjc,dbesyc) 
else 

use asymptotic approximation 

call cbesjy 3(arg, m, icode, besjc,besyc,dbesjc,dbesyc) 
endif 




return 

end 

subroutine chankel (arg, order, icode, hanic,han2c,dhanlc,dhan2c) 

c This set of routines calculates various types of Hankel functions 
c in single precision, 

c INPUTS: 


c 

c 

c 

c 

c 


arg — complex argument 

order — integer order m >= 0 

icode -- 1 calculate H[l]m(x), the Hankel fen of first kind 

2 Il[2]m(x), the Hankel fen of second kind 

3 H[l]m(x),H[2]m(x) & first derivatives 


c OUTPUT: 


c the value of the indicated Hankel function, written to the following 
c variables, according to the value of icode: 
c icode = I — > hanlc 

c 2 han2c 

c 3 hanlc, han2c,dhanlc,dhan2c 


integer order, icode 

complex arg, hanlc, dhan Ic,han2c,dhan2c 
complex besjc,dbesjc,besyc,dbesyc,ic 

ic = (0.,1.) 

call cbessel(arg,order,3,besjc,besyc,dbesjc,dbesyc) 
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if (icode .eq. 1) then 

haulc = besjc + ic*besyc 
else if (icode .eq. 2) then 
han2c = besjc - ic*besyc 
else 

hanlc = besjc + ic*besyc 
han2c = besjc - ic*besyc 
dhanlc = dbesjc + ic*dbesyc 
dhan2c = dbesjc - ic*dbesyc 
endif 

return 

end 

******************************** ********** 
subroutine cbsjl(arg,m,besjl) 

integer m 
complex arg,besjl 
double complex darg,besj 
double precision ertolp,ertolm 

i 

common /DERRTOL/ ertolp,ertolm 
darg — arg 

call dcbsji(darg,in,besj) 
besjl = besj 

return 

end 

^******************************************* 
subroutine cdbsjl(arg,m,dbesjl) 

integer m 

complex arg,dbesjl 
double complex darg,dbesj 
double precision ertolp,ertolm 

common /DERRTOL/ ertolp,ertolrn 

darg — arg 

call dcdbsjl(darg,m,dbesj) 
dbesjl = dbesj 

return 

end 
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^******%********************* *************** 
subroutine cbsy 1 (arg, rn, besy 1) 

integer m 
complex arg, besy 1 
double complex darg,besy 
double precision ertolp,ertolm 

common /DERRTOL/ ertolp,erto!m 

darg — arg 

call dcbsyl(darg,m,besy) 
besyl = besy 

return 

end 

c * **********+**+***+************++*********+ 
subroutine cdbsyl(arg,m, dbesy 1) 

integer m 

complex arg, dbesy 1 
double complex darg, dbesy 
double precision erto!p,ertolm 

common /DERRTOL/ ertolp,ertolm 

darg “ arg 

call dcdbsyl (darg, m, dbesy) 
dbesyl == dbesy 

return 

end 

subroutine cbesjy3(arg,m,icode,besj3,besy3,dbesj3,dbesy3) 
c Airy function approximation to Bessel function, complex argument 
integer m,icode 

complex arg,besj3,besy3,dbesj3,dbesy3 
double complex darg, besj, besy ,dbesj, dbesy 
double precision ertolp,ertolm 

common /DERRTOL/ ertolp,ertolm 

darg = arg 

call dcbesjy3(darg, m,icode, besj, besy, dbesj, dbesy) 
besj 3 = besj 
besy3 = besy 
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dbesj 3 = dhesj 
dbesy3 = dbesy 

return 

end 

c **********+******** ***+***+******#********* 

subroutine cbesjy2(arg,m,icode,besj2,besy2,dbesj2,dbesy2) 

c asymptotic approximation to Bessel function, complex argument 
c Abr-Stegun 9.3.35, 9.3.36 

integer m,icode 

complex arg, besj 2, besy 2, dbesj2, dbesy 2 
double complex darg,besj,besy,dbesj, dbesy 
double precision ertolp,ertolm 

common /DERRTOL/ ertolp,ertolm 

darg = arg 

cal 1 dcbesjy2(darg,m, icode, besj , besy , dbesj , dbesy ) 
besj2 = besj 
besy2 = besy 
d besj 2 = dbesj 
dbesy2 = dbesy 

return 

end 

subroutine rdbessel(arg,m,icode, besj, besy, dbesj, dbesy) 

c This set of routines calculates various types of Bessel functions 

c in double precision. 

c INPUTS: : - 

c arg -- real argument x >= 0 
cm— integer order m >— 0 

c icode — 1 calculate Jrn(x), the Bessel fen of first kind 
c 2 Ym(x), the Bessel fen of second kind 

c 3 Jm(x), Ym(x) and their first derivatives 

c OUTPUT: 

c the value of the indicated Bessel function, written to the following 
c variables, according to the value of icode: 
c icode = 1 — > besj 

c 2 besy 

c 3 besj, besy, dbesj, dbesy 
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c NOTE: Yrn(O) and Y’m(O) are undefined. 


c Three different evaluation methods are used, depending on the 
c relationship between M and ARG: 

c 1) if both are small, the power series definition is used directly 
c (routines DRBSJl, DRDBSJl, DRBSYl, DRDBSYl). 

c 2) if both are large, a series involving Airy functions is used 
c (routine DRBESJY3). 

c 3) if M is small but ARG is large, an asymptotic expansion is used 
c (routine DRBESJY2). 

integer m,icode 

double precision arg,zm,besj,dbesj,besy,dbesy,ertolp,ertolm 
double precision zarg 

common /DERRTOL/ ertolp,ertolm 

ertolp ~ 1.0d35 
ertolm = 1.0d-35 

zarg = min( .007*m**2+.16*m+12.5, 40. ) 
zm = min( 3.*arg/ll.+20., 33. ) 

if (arg .le. zarg) then 
c use power series 

if (icode .eq, I) call drbsjl(arg,rn,besj) 
if (icode .eq. 2) call drbsy l(arg,m,besy) 
if (icode .eq. 3) then 
call drbsjl(arg,m,besj) 
call drdbsjl(arg,m,dbesj) 
call drbsy 1 (arg, m,besy) 
call drdbsy l(arg,m,dbesy) 
endif 

else if (m .ge. zm) then 
c use airy function approximation 

call drbesjy2 (arg, m, icode, besj,besy,dbesj,dbesy) 
else 

c use asymptotic expansion approximation 

call drbesjy3(arg,m,icode,besj,besy,dbesj,dbesy) 
endif 

return 

end 

subroutine rdhankel(arg,order,icode,hanl,han2,dhanl,dhan2) 

c This set of routines calculates various types of Hankel functions 
c in double precision. 
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c INPUTS: 


c arg — real argument x >= 0 
c order — integer order m >= 0 

c icode — 1 calculate H[l]m(x), the Hankel fen of first kind 

c 2 H[2]in(x), the Hankel fen of second kind 

c 3 H[l]m(x), II[2]m(x) & first derivatives 

c OUTPUT: 

c the value of the indicated Hankel function, written to the following 
c variables, according to the value of icode: 
c icode = I --> hanl 
c 2 han2 

c 3 hanl,han2,dhanl,dhan2 

c NOTE: the variables hanl,han2,dhanl,dhan2 must be declared complex, 

integer order, icode 

double complex hanl,dhanl,han2,dhan2,ic 
double precision arg, besj, dbesj, besy,dbesy 

ic = (0.d0,l.d0) 

if (arg .ge. O.dO) then 

call rdbessel(arg,order,3,besj,besy,dbesj,dbesy) 
if (icode .eq. 1) then 
haul = besj + ic*besy 
else if (icode ,eq. 2) then 
han2 = besj - ic*besy 
else 

hanl = besj + ic*besy 
han2 = besj - ic*besy 
dhanl — dbesj + ic*dbesy 
dh an 2 = dbesj - ic*dbesy 
endif 
else 

arg = -arg 

c Abr-Stegun 9.1.39 

call rdbessel (arg, order, 3, besj, besy, dbesj ,dbesy) 
if (icode .eq. 1) then 

hanl = -(-l.dO)**order*(besj-ic*besy) 
else if (icode .eq. 2) then 

han2 = -(-l,dO)**order*(besj fic*besy) 
else 

hanl = -(-l.dO)**order*(besj-ic*besy) 
hari2 = -(-l.doj^^order^lbesj+ic^besy) 
dhanl = (-l.dO) + 4: order*(dbesj-ic*dbesy) 
dhan2 = (-l.dO)**order*(dbesj+ic*dbesy) 
endif 
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endif 


return 

end 


c * ****************+************************* 
subroutine drbsjl(arg,m, besjl) 

integer rn,lone,k,mone 

double precision arg,besj 1 ,fm,qntone,fmone,f,rdfact, ertolp, ertolm 

common /DERRTOL/ ertolp, ertolm 

f = rdfacl(m,arg) 
if (f .It. 1.5d0*ertolin) then 
besjl = O.OdO 
return 
endif 

if (f ,gt. .5d0*ertolp) then 
besjl = ertolp 
return 
endif 

fm = dble(m) 

lone = int(10.0d0 + 4.0d0*arg/3.0d0) 
qntone = l.OdO 

do 10 k = l, lone 
inone = lone - k + 1 
fmone = dble(mone) 

qntone = l.OdO - qntone*(0.5d0*arg)**2/(fmone*(fmone + fm)) 
10 continue 

besjl = f*qntone 

return 

end 

function rdfact(n,arg) 
c calculate (arg/2) A n/n! 
integer n 

double precision rdfact,arg, flag, ertolp, ertolm 

common /DERRTOL/ ertolp, ertolm 

rdfact = l.OdO 
if (n .eq. 0) return 
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do 5 k = l,n 

rdfact = rdfact*arg/(2.0d0*dble(n~k+l)) 
if (rdfact .it. ertolm) then 
rdfact = ertolm 
return 
endif 

if (rdfact .gt. ertolp) then 
rdfact = ertolp 
return 

endif .... . 

5 continue 
return 
end 

subroutine drdbsjl(arg,m, dbesjl) 
integer m, lexit, k,mexit 

double precision arg, dbesjl, frn,dqnton,amexit,f, rdfact 
double precision ertolp, ertolm, besjl 

common /DERRTOL/ ertolp, ertolm 

if (m .eq. 0) then 

call drbsj 1 (arg, 1, besjl) 
dbesjl = -besjl 
return 
endif 

f = rdfact (m,arg) 
if (f .It. 1.5d0*ertolm) then 
dbesjl = O.OdO 
return 
endif 

if (f .gt. .5d0*ertolp) then 
dbesjl = ertolp 
return 
endif 

fm — dble(m) 

lexit = int(i0.0d0 + 1.6d0*arg) 
dqnton = l.OdO 
do 10 k = 1, lexit 

mexit = lexit - k + 1 
amexit = dble(rnexit) 

dqnton = l.OdO - dqnton*(0.5d0*arg)**2*(fm+2.0d0*amexit) / 
. (amexit* (frn+amexit)*(frn+2.0d0*arnexit-2.0d0)} 

10 continue 

dbesjl = 0.5d0*dqnton*rdfact(m-l,arg) 
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return 

end 


c * ****************************************** 
subroutine drbsyl(arg,m,besyl) 

integer m,k,lz,klz,iflag 

double precision arg,besyl,fm, gamma, pi, repi,emz,besjl,pemzj, qntbz 
double precision am,sumdk,gmz 

double precision rdfact,f,fmz,ertolp,ertolm,qntdk,fklz,dk,ddigam 

common /DERRTOL/ ertolp,ertolm 

f zz: rdfact(m,arg) 
fm = dble(m) 

gamma = 0.577215664901 53d0 
pi = 3.1415926535897932d0 
repi — l.OdO/pi 

call rdtolch (arg,besy 1 ,0.0d0,-ertolp,iflag) 
if (iflag .eq. 1) return 

emz = 2.0d0*repi*(gamma + dlog(0.5d0*arg)) 
call drbsjl(arg,m,besjl) 
pemzj — emz*besjl 

if (in .eq. 0) then 
fmz — O.OdO 
else if (m .eq. 1) then 
fmz = -2.d0*repi/arg 
else if (f .It. 1.5d0*ertolm) then 
fmz = -ertolp 

else if (f .gt. .5d0*ertolp) then 
fmz = O.OdO 
else 

qntbz = l.OdO 
do 10 k = 1 ,m-l 

am = dble(m-k) 5 . > ' ■ 

qntbz = l.OdO + qntbz*(0.5d0*arg)**2/(am*(fm-am)) 

10 continue 

fmz = -qntbz*repi / (f*fm) 
endif 

lz = inl(7.0d0+1.5d0*arg) 
qntdk = l.OdO 
do 20 k = 1,1® 
klz = lz - k + 1 
fklz = dble(klz) 

dk = (ddigam(klz+l) + ddigam(m+klz+l)) / 

((ddigain(klz)+ddigam(klz+m))*(fklz+l.)*(fm+fklz+l.)) 
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qntdk = l.OdO - qntdk*dk*(0.5d0*arg)**2 
20 continue 

sumdk = ddigam(m)-(qntdk*(l.d0+ddigam(m+l))*(arg*0.5d0)**2 / 
(fm+1.0d0)) 

if (f .It. 1.5dO*ertolm) then 
gmz = O.OdO 

else if (f .gt. ,5d0*ertolp) then 
gmz — ertolp 
else 

gmz = -sumdk*repi*f 
endif 

besyl = pemzj + fmz 4- gmz 

return 

end 

$ $ $ $ $ ^ 4 : ^ ^ $ $ $ $ $ $ $ £ £ £ * $ £ - 1 ^ $ $ $ $ $ $ $ :£ 

function ddigam(n) 
integer n j 

double precision ddigam 

ddigam = O.OdO 
if (n ,eq. 0) return 
do 5 j = n,l,-l 

ddigam = ddigam + l.d0/dble(j) 

5 continue 

return 

end 

^$ $ $ $ $ $ 4 $ $ 4 $ 4 4 $ 4 4 : 4 ’ 4 $ $ 4 $ $ 4 4 $ 4 $ $ $ 4 $ $ 4 $ $ $ $ $ 4 $ 

subroutine drdbsy l(arg,m,dbesyl) 
integer m,k, mm, kmm,larg,klarg, iflag 

double precision arg,dbesyl,fm, gamma, pi, repi,emz,dbesjl,dcmzjl 
double precision dernzj2,qntdb,f,fk,dfmz,qntder,flarg,dkder,ddigain 
double precision dergmz,ertolp,ertolm,rdfact,besjl,sumder 

common /DERRTOL/ ertolp, ertolm 

f = rdfact(m,arg) 
fin = dble(rn) 

gamma — 0.57721566490153d0 
pi = 3.1415926535897932d0 
repi = l.OdO/pi 

call rdtolch(arg,dbesy 1, O.OdO, ertolp, iflag) 
if (iflag .eq. 1) return . . . 
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emz = 2.0d0*repi* (gamma + dlog(0.5d0*arg)) 

call drdbsji (arg,m,dbesjl) 

demzjl = emz*dbesjl 

call drbsjl(arg,m,besjl) 

deinzj2 = besji*2,OdO*repi/arg 

if (in .eq. 0) then 
dfmz = O.OdO 
else if (m .eq. 1) then 

dfmz = 2.0d0*repi/arg**2 
else 

qntdb = l./(.5d0*arg*f) 
k = 1 

mm = m - 1 
10 kinm = mm - k 
fk = dble(k) 

qntdb = qntdb + ((fm-(2.dO t fk)) + rdfact(k,arg)) / 
(rdfact(kmm,arg)*(.5d0*arg)**2) 
if (abs(qntdb) .It. 1.5d0*ertolrn) then 
dfmz = O.OdO 

else if (qntdb .gt. .5d0*ertolp) then 
dfmz — ertolp 
else 

dfmz = .5d0*qntdb*repi 
end if 
k = k + 1 

if (k .le. mm) go to 10 
endif 

if (f .It. 1.5d0*ertolm) then r - ; 

dergmz — O.OdO 
else if (f .gt. .5d0*ertolp) then 
dergmz = ertolp 
else 

larg = int(8.0d0+1.4d0*arg) 
qntder = l.OdO 
do 20 k — l,larg 

klarg = larg - k + l 
flarg = dble(larg~k4d) 

dkder = (ddigam(klarg-f-l) + ddigam(m+klarg+l)) * 
(fm4-2.d0*flarg+2.d0) 

/ ( (ddigam(klarg)+ddigam(nH*klarg)) + (flarg -fl.dO)* 
(fm+flarg + I.d0)*(fm+2.d0*flarg) ) 
qntder = l.OdO - qntder*dkder*(0.5d0*arg)**2 
20 continue 

sunider = fm*ddigam(m) - qntder*(l.dO+ddigam(m+ 1))* 

(fm 4- 2.d0) + (0.5d0*arg)**2/(fm-f l.dO) 
dergmz = -0.5d0*repi*f*sumder/(0.5d0*arg) 
endif 
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dbesyl = demzjl + demzj2 + dfmz + dergmz 
call rdtolch(dbesyl,dbesyl ? ertolp,O.OdO,iflag) 
return 
end 

subroutine drbesjy2(arg,m,icode,besj2,besy2,dbesj2,dbesy2) 
c Airy function approximation to Bessel function, real argument 
integer rn,icode 

double precision arg,besj2,besy2,dbesj2,dbesy2,fm,x,zeta,ql,phi 
double precision qli,g3,g32,g2,gsq,al,bO,cO,dl,h,y,airy,biry 
double precision ertolp,ertolm,fl,f2,f3, dairy ,dbiry 
double precision ul,u2,u3,bl,q,u4,u5,a2,b2 

common /DERRTOL/ ertolp,ertolm 

fm = dble(m) 

x = arg/fm 

if (x .gt. l.OdO) then 

zeta = -(1.5d0*(sqrt(x**2-l.d0) - acos(l.d0/x)))**(2.d0/3.d0) 
else 

zeta = (I.5d0*(dlog((1.0d0 +sqrt(l.d0-x**2))/x) 

- sqrt(l.d0-x**2)))**(2.d0/3.d0) 

endif 

if (x .gt. 0.98d0 .and. x .It. 1.02d0) then 
h = l.OdO - x 

phi = 2.0d0**(l.d0/3.d0)*(l.0d0 + 0.2d0*h + 3.d0*h**2/35.d0 
+ 73.dO + h**3/1575.dO + 35209.d0*h**4/l212750.d0 + 
380069.d0*h + *5/l8768750.d0) 
al = -1. dO/225. dO - 7i.d0*h/38500.d0 
bO = (l. d0/70. dO + 2.dO + h/225.dO)*2.dO**(l.dO/3.dO) 
bl = 2.dO**(l.dO/3.)*(-1213.dO/1023750.dO - 3757.d0*h/2695000. 
- h**2*(8.9962899979797d-4 f h*(.0002753433716d0 - h* 
(.00018048868d0 + h*.0004 108523)))) 
a2 = 6.93735541 3546877d-4 + h*(.00046448349036601 - h* 
(.0002890362546053d0 + h*(.0008747649439535d0 + 
h*. 0010297 1637614))) 

b2 = -2.d0**(l./3.)*(4.382918094497229d-4 + h* 

(7.1 1048651 1691 ld-4 + h*5.318984348085d-4)) 
b2 = b2 - 2 .dO* * ( 1 . / 3 . ) * h * * 3 * 2 . 1 8 295 8472d- 4 
cO = (O.ldO + 0.02d0*h)*2.d0**(2.d0/3.d0) 
dl = 23. d0/3150. dO -f 1453.d0*h/346500.d0 

else 

ql = zeta/(l.0d0-x**2) 
phi = (4.0d0*ql)**0.25d0 
fl = l.OdO - x**2 
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12 = fl**2 ' ’ - 1 

f3 = fl**3 

qli = l.OdO/ql 

g3 = qti**3 

g32 = qli**1.5dO 

g2 = qli**2 

gsq = sqrt (qli) 

al = ( -455.d0*g3/4608.d0-7.d0*g32*(fl-5.d0/3.d0)/384.d0 + 
385.dO/1152.dU - 77.dO*fl/J92.dO + 9.dO*f2/128.dO ) / f3 
bO = (-5.dO*g2/48.dO + gsq*(5.dO/24.dO - fl/8.d0))/f2 
cO = (7.dO*qli/48.dO - 3.dO*sqrt(ql)*(7.dO/72.dO - fl/8.dO))/fl 
dl = ( 385.d0*g3/4608.d0 + 5.dO*g32*(-rt + 7.dO/9.dO)/128.dO - 
15.dO*f2/128.dO + 33.d0 *fl/64.d0 - 455.dO/1152.dO ) / 13 


q = sqrt,(fl/zeta) 

ul = q*(i.dO-5.dO/(fl*3.dO)) / (8.dO*fl) 

u2 = (4.5dO - 77.dO/(fl*3.dO)*(l.dO-5.dO/(fl*6.dO))) / 

(64.d0*fl) 

u3 = q*(75.dO/2. - 456.3dO/fl + 17017.d0/(f2*18.d0)*(l.d0- 
5.dO/(fl*9.dO))) / (512.dO*f2) 

u4 = (3675.dO/8.dO - 9683.3d0/fl + 2717.dO/f2*(53.dO/4.dO - 
2737.dO/(fl *162.d0)*( 1 ,d0-5.d0/(12.d0*fl))))/(4096.d0*f2) 
u5 = 59535. dO/8.dO - 221.d0/(4.d0*fl)*(305923.d0/70.d0 - 
77. dO/fl *(14743. dO/45. dO - 95.dO/fl*(67.dO/9.dO- 
3335.d0/(486.d0*fl)*(l.d0-l.d0/(3.d0*fl))))) 
u5 = q*u5/(f3*32768.d0) 

bl = -u3 - 5.d0/(zeta**2*48.d0)*(u2+77.d0/(zeta*96.d0)* 

(ul + 221.dO/(zcta**2*144.dO))) 
b2 = -u5 - 5.d0/(zeta**2*48.d0)*(u4 + 77.d0/(zeta*96.d0)* 

(u3 + 221.dO/(zeta**2*144.dO)*(u2 + 437.dO/(zeta*192.dO)* 
(ul + 145. dO/ (zela**2*48.dO))))) 
a2 = u4 - 7.dO/(zeta*48.dO)*(u3 + 65.d0/(zeta**2*96.d0)*(u2 + 
209.d0/(zeta*144.d0)*(ul + 425.dO/(zeta**2*192.dO)))) 

endif 


y = zeta*fm**(2.dO/3.dO) 

call rdairy fn(y,icode, airy, biry, dairy ,dbiry) 


if (icode .eq. 1) ( hen 

besj2 = phi/(fm**(l.dO/3.dO))*(airy*(l.dO+al/fm**2-ta2/fm**4) 
+ dairy / (fm**(4.dO/3.dO))*(bO+bl/fm**2+b2/fm**4)) 
else if (icode .eq. 2) then 

if (biry .eq. ertolp .or. dbiry .eq. ertolp) then 
besy2 = -ertolp 
else 

besy2 = -phi/(fm**(l./3.))*(biry*(l.d0+al/fm**2+a2/frn**4) 
+ dbiry / (fm**(4.dO/3.dO))*(bO-t-bl/fm**2+b2/fm**4)) 

endif 

else 
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besj2 - phi/(fm^(l.dO/3.dO))*(airy*(l.dO+al/fm**2+a2/fni**4) 
+ dairy / (fm**(4.dO/3.dO))*(bO-hbl/fm**2 f b2/fm**4)) 


if (biry .eq. ertolp .or. dbiry .eq. ertolp) then 
hesy2 = -ertolp 
else 

besy2 — -phi/(fm**(l./3.))*(biry*(l.d0+al/fm**2 fa2/fm**4) 
+ dbiry / (fm**(4.d0/3.d0))*(b0+bl/fin* :,t 2 h b2/fm**4)) 

endif 

dbesj2 = -2,d0/( rm + *(2.dO/3.dO) + x + phi)*(airy/(fm**(2.dO/3.dO)) 
*c0 + dairy*(i.dO + dl/(fm**2.d0)) ) 

if (biry .eq. ertolp .or. dbiry .eq. ertolp) then 
dbesy2 = ertolp 
else 

dbesy2 = 2.dO/(rni**(2.dO/3.dO)*x*ph0*(Wry*cO/ 

(fm**(2.dO/3.dO)) + dbiry*(l.dO + dl/(fm**2.d0))) 
endif » 

endif 

return 

end 

subroutine rdairy fn(y,icode, airy, biry , dairy, dbiry) 
integer icode,ll,k,k2,l3,l4,l5,m5,Ip,lpt 

double precision y, airy, biry, dairy, dbiry, cl, c2,flam,reflam, pi 
double precision dam,srtwo,qntv,ffml,vffml,vffmr,qntw,ffm2,wffmr 
double precision qntdv,ffm3,dvffmr,qntdw,ffm4,dwffmr,surni,sum2 
double precision sum4,diftl,sumtl,dift2,surnt2,qntf,qnfcg,qntdf 
double precision ff,fg,fdf,fdg,frnl,fpl,fmp,gml,gpl,frnpt,sxfmpt 
double precision ertolp, ertolm,srpi,aby , slam, sum3,ffrn5,qntdg 
double precision sxfmp 

common /DERRTOL/ ertolp, ertolm 

cl = 0, 3550280538878 17d0 
c2 = 0.2588 19403792807d0 

if (y .gt. -lO.OdO .and. y .It. 6.0d0) then 

c power series 

if (y *6 e * O.OdO) then 
L5 = int(5.0d0 4* 3.0d0*dabs(y)) 
else 

L5 = int(5.0d0 + 3.5d0*dabs(y)) 
endif 

qntf = l.OdO 
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qntg = 1.0(10 
qntdf = l.OdO 
qntdg = l.OdO 
do 84 k = 1,L5 
m5 = L5 - k + 1 
ffm5 = 3.dO*dble(m5) 

qntf = l.OdO + qntf*y**3/(ffm5*(ffm5 - l.OdO)) 
qntg = l.OdO + qntg*y**3/(ffm5*(ffm5 + l.OdO)) 
qntdf — l.OdO + qntdf*y**3/(ffni5*(ffm5 + 2.0d0)) 
qntdg = l.OdO + qntdg*y**3/(ffm5*(ffm5-2.0d0)) 
84 continue 
ff = qntf 
fg = qntg*y 

fdf = qntdf*0.5d0*y**2 

fdg = qntdg 

airy = cl*ff - c2*fg 

dairy = cl*fdf - c2*fdg 

biry = dsqrt(3.0d0)*(cl*ff + c2*fg) 

dbiry = dsqrt(3.0d0) + (cl*fdf + c2*fdg) 


flam — (2.dO/3.dO)*dabs(y)**i.5dO 

ref] am = l.OdO/flam 

pi = 3.1415926535897932dO 

srpi = dsqrt(pi) 

aby = dabs(y)**0.25d0 

if (y .gt. 25. dO) then 
airy = O.OdO 
dairy = O.OdO 

biry = ertolp ^ *■ •• 

dbiry — ertolp 

else if (y .ge. 6.0d0 .and. y .le. 25.0d0) then 

c exponential approximation 

lp = int(4.d0 + 300.0d0*reflam) 

if (reflam .gt. .06) lp = 22 + ll*min((reflam-.06)/.003, 

(.102-reflarn) /.039) 

fml = l.OdO 
fpl = l.OdO 
do 115 k = l,lp 

sxfmp = 6.0d0*dble(lp-k+l) 

fml = l.dO - fnil*reflani*(sxfmp-5.d0)*(sxfmp-l.d0) 
/(12.d0*sxfmp) 

fpl = l.dO j- fpl*reflam*(sxfmp-5.dO)*(sxfmp-l.dO) 
/(12.d0*sxfmp) 

115 continue 
Ipt = lp 
gml — l.OdO 
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gpl = l.OdO 
do 125 k = ljpfc 

sxfmpt = 6.0d0*dble(lpt-k+l) 

gml = LdO-gml*reflam*(sxfmpt+l.dO)*(sxfmpt-7,dO) 
/(22,d0*sxfmpt) 

gpl — I.d0+gpl*reflam*(sxfmpt+l.d0)*(sxfmpt-7.d0) 
/(12.d0*sxfmpt) 

125 continue 

airy = 0.5d0*y t *(“0.25d0)*fml*exp(-flam)/9rpi 
dairy — -0.5d0*y**0.25d0*gml*exp(-flam)/srpi 
biry = y**(d).25dO)*fpl*exp(flam)/srpi 
dbiry = y**0.25d0*gpl*exp(flam)/srpi 

else 

c trigonometric approximation 

slam = sin(flarn) 
clam = cos(flarn) 
srtwo = sqrt(2.d0)/2.d0 
Ll = irit (2-0d0 + 160.0d0*reflam) 
qntv = l.OdO 
do 20 k = 1,L1 
ffml = 12.d0*dble(Ll-k+l) 

vffinr — (ffrnl-1 l.dO)*(ffnil-7.dO)*(fTinl-5.dO)*(ffml-l.dO) 

/ (I44.d0*ffml*(ffinl - 6.d0)) 
qntv = l.OdO - (qntv*vffmr*reflam**2.d0) 

20 continue 

L2 = int(2.0d0 + 220.0d0*reflam) 
qntw = l.OdO 
do 30 k = 1,L2 
ffm2 = 12.dO*dble(L2-k+l) 

wffmr = (ffm2-5.dO)*(ffm2-l.dO)*(frm2+l.dO)*(ffm2+5.dO) 

/ (I44.d0*ffm2*(ffm2 + 6.dO)J 
qntw = l.OdO - (qntw*wffmr*reflam**2.d0) 

30 continue 

L3 = int(2.0d0 + 130.0d0*reflam) 
qntdv = l.OdO 
do 40 k = l,L3 
ffm3 = 12.d0*dble(L3-k+l) 

dvffmr - (ffm3-7.d0)*(ffm3-l.d0) + (ffm3+Ld0)*(ffm3+7.d0) 
/ ( 1 44 . dO * ffm3 * ( f fm 3 + 6.d0)) 
qntdv = l.OdO - (qntdv *dvffmr*reflam**2.d0) 

40 continue 

L4 = int(2.0d0 + 200.0d0*reflam) 
qntdw = l.OdO 
do 50 k = 1,L4 
ffrn4 = 12,dO*dble(L4-k+l) 

dwffmr = (ffm4-13,d0)*({fm4-7.d0)*(ffm4-5.d0) t (ffm4+Ld0) 
/ (I44,d0*ffin4*(ffm4 - 6.d0)) 
qntdw — l.OdO - (qntdw*dwffmr*reflam**2.d0) 
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50 continue 

surnl = qntv/srpi 

sum2 = 5.d0*qntw*reflam/(72.d0*srpi) 

suniS — -72.d0*qntdv*reflam/(7.d0*srpi) 

sum4 = qntdw/srpi 

diftl = srtwo*(suml - sum2) 

sumtl = srtwo*(surnl + sum2) 

dift2 = srtwo*(sum4 - sum3) 

sumt2 — srtwo*(sum4+sum3) 

airy — (diftl*clam + sumtl*slam)/aby 
dairy = (dift2*slam - surnt2*darn)*aby 
biry = (sumtl *clam - diftl + slam) /aby 
dbiry = (sumt2*slam + dift2*clam) i aby 

endif 

cndif 

return 

end 

c ^^ + ^5f: + + * + * + **** + **** + ****************** : t : i : ** 

subroutine drbesjy3(arg,m,icode,besj3,besy3,dbesj3,dbesy3) 
c asymptotic approximation to Bessel function, real argument 
integer m,icode,lmk,k 

double precision arg,besj3,besy3,dbesj3,dbesy3,fm,pi,srarg,angle 
double precision qnta,qntb,f]mk4,amk,bmk,umk,vrnk,qntad,qntbd,arnkd 
double precision smk,ertolp,ertohn,fmsr,bmkd,rmk 

common /DERRTOL/ ertolp,ertolm 

frn = dble(m) 

pi = 3.1415926535897932d0 

srarg = sqrt(2,0d0/(pi*arg)) 

angle = arg - 0.25d0*pi*(1.0d0 + 2.0d0*fm) 

finsr = 4.0d0*fm**2 

lmk = int( 2 -| m/3. + 27./sqrt(arg-12.) ) 
if (link .It. 1) link = 1 

qnta = l.OdO 
qritb = J .OdU 
do 10 k = l,lmk 

flrnk4 = 4.0d0*dble(lmk-k-| 1 ) 

amk — (4.0d0*(fmsi-(flink4-3.0d0)**2)*(fmsr-(flmk4-1.0d0)**2)) / 
(flmk4*(flmk4 - 2.0d0)) 

bmk = (4.0d0*(fmsr-(flmk4-1.0d0)**2)*(fmsr-(flmk4+1.0d0)**2)) / 
(flmk4*(flmk4 + 2.0d0)) 
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qnta — l.OdO - qnta*amk/(8.0d0*arg)**2 
qnfcb = l.OdO - qntb*bmk/(8.0d0*arg)**2 
10 continue 
umk — qnta 

vmk — qntb*(fmsr - 1.0d0)/(8.0d0*arg) 
besj3 = srarg*(umk*cos(angle) - vmk*sin(angle)) 
besy3 = srarg*(umk*sin(angle) + vmk*cos(angle)) 
call rdtolch(besj3,besj3,ertolp,0.0d0, iflag) 
call rdtolch(besy3,besy3,ertolp,0.0d0,iflag) 

if (icode .eq. 3) then 
qntad = l.OdO 
qntbd = l.OdO 
do 20 k = l,lmk 
flmk4 = 4.0d0*dble(lmk-k+l) 
amkd = (4.0d0*(fmsr+(nmk4**2-1.0d0))* 

(fmsr-(nmk4-3.0d0)* + 2) + (fmflr-(flmk4-5.0d0) ? * :t 2)) / 
(flmk4*(flmk4 - 2.0d0)*(fmsr-f (flmk4-4.0d0)**2-1.0d0)) 
bmkd = (4.0d0*(fmsr+(flmk442.0d0)**2-1.0d0)* 

(fmsr-(llnik4-3.0d0)**2)*(fmsr-(flmk4-1.0d0)**2)) / 
(flmk4*(flmk4+2.0d0)*(fmsr+(flmk4-2*0d0)**2-1.0d0)) 
qntad = l.OdO - qntad*amkd/(8.0d0*arg)**2 
qntbd = l.OdO - qntbd*bmkd/(8.0d0*arg)**2 
20 continue 

rink = qntad 

sink = qntbd* (fmsr 4~ 3.0d0)/(8.0d0*arg) 
dbesj3 = -srarg*(rmk*sin(angle) + smk*cos(angle)) 
dbesy3 = srarg*(rmk*cos(angle) - snik*sin(angle)) 
call rdtolch(dbesj3,dbesj3,ertolp,O.OdO, iflag) 
call rdtolch(dbesy3,dbesy3,ertolp,0.0d0, iflag) 
endif 

return 

end 


£**************+**************************** 
subroutine rdtoIch(rin, rout, setbig,selsmall, iflag) 

c real, double precision check to avoid over- and underflow 

double precision rin,rout,setbig,setsmall,ertolp,ertolm 
integer iflag 

common /DERRTOL/ ertolp,ertolm 
iflag = 0 

if (dabs(rin) .It. 1.5d0*ertolm) then 
rout = setsmall 
iflag = 1 

endif , * ^ 

if (dabs(rin) .gt. .5d0*ertolp) then 
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rout = setbig 
iflag = 1 
endif 
return 
end 

*****************************************+ 

subroutine cdbessel (arg, m, icode, besjc, besyc,dbesjc,dbesyc) 

c This set of routines calculates various types of Bessel functions 
c in double precision. 

c INPUTS: 

c arg — complex argument 
c m — integer order m >= 0 

c icode — 1 calculate Jrn(x), the Bessel fen of first kind 

c 2 Ym(x), the Bessel fen of second kind 

c 3 Jm(x),Ym(x) & their first derivatives 

c OUTPUT: 

c the value of the indicated Bessel function, written to the following 
c variables, according to the value of icode: 
c icode = 1 --> besjc 
c 2 besyc 

c 3 besjc, besyc, dbesjc,dbesyc 

c NOTE: Ym(0) and Yhn(O) are undefined. 

c Three different evaluation methods are used, depending on the 
c relationship between M and ARG: 

c 1) if both are small, the power series definition is used directly 
c (routines DCBSJi, DCDBSJl, DCBSY1, DCDBSYI). 

c 2) if both are large, a series involving Airy functions is used 
c (routine DCBESJY3). 

c 3) if M is small but ARG is large, an asymptotic expansion is used 
c (routine DCBESJY2). 

integer m, icode 

double complex arg, besjc, besyc, dbesjc,dbesyc 
double precision ertolp,ertolm,zm 

common /DERRTOL/ ertolp,ertolm 

ertolp = 1.0d35 
ertolrn = I.0d-35 

if (m .It. 8) then 

zm = 7. dO H- 2.d0*dble(m)/7.d0 
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else 

zm = (184, dO - 8.dO*dble(rn))/15.dO 
endif 

if (zabs(arg) .le. zm) then 
c use power series 

if (icode .eq. l) call dcbsjl(arg,m,besjc) 
if (icode .eq. 2) call dcbsyl(arg,m,besyc) 
if (icode .eq. 3) then 
call dcbsjl(arg,m,besjc) 
call dcdbsjl(arg,m,dbesjc) 
call dcbsyl (arg,m,besyc) 
call dcdbsyl(arg,m,dbesyc) 
endif 

else if (m .ge. 8) then 
c use airy approximation 

call dcbesjy2(arg,m, icode, besjc,besyc,dbesjc,dbesyc) 
else 

c use asymptotic approximation 

call dcbesjy3(arg,m, icode, besjc,besyc,dbesjc,dbesyc) 
endif 

return 

end 

subroutine cdhankel(arg,order,icode,hanlc,han2c,dhanlc,dhan2c) 

c This set of routines calculates various types of Hankel functions 
c in double precision. 

c INPUTS: 

c arg — complex argument 
c order — integer order m >= 0 

c icode — 1 calculate II[l]in(x), the Hankel fen of first kind 
c 2 H(2]m(x), the Hankel fen of second kind 

c 3 H[l]in(x),H[2]in(x) & first derivatives 

c OUTPUT: 

c the value of the indicated Hankel function, written to the following 
c variables, according to the value of icode: 
c icode = 1 --> hanlc 

c 2 han2c 

c 3 hanlc, han2c,dhanlc,dhan2c 


integer order, icode 

double complex arg, hanlc, dhanlc,han2c,dhan2c 
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double complex besjc, dbesjc, besyc,dbesyc,ic 
ic = (0.d0,l.d0) 

call cdbessel(arg, order, 3, besjc, besyc, dbesjc, dbesyc) 
if (icode .eq. 1) then 

hanlc = besjc + ic*besyc 
else if (icode .eq. 2) then 
han2c = besjc - ic*besyc 
else 

hanlc = besjc + ic*besyc 
han2c = besjc - ic* besyc 
dhanlc = dbesjc + ic*dbesyc 
dhan2c — dbesjc - ic*dbesyc 
endif 

return 

end 


^* *********************+************ ******** 
subroutine dcbsjl(arg,m, besjl) 

integer in, lone, k 

double complex a.rg,besjl,qntone,f,cdfact,ic 
double precision ertolp,ertolm,fm,fmone, theta 

common /DERRTOL/ ertolp,ertolm 

f = cdfact(m,arg) 

fm — dble(m) 
ic = (0.dU,l.d0) 

if (zabs(f) .It. 1.5d0*ertolm) then 
besjl — O.OdO 

else if (zabs(f) .gt. .5d0*ertolp) then 
theta = datan(dimag(f)/dble(f)) 
besjl — ertolp*exp(ic*fm*theta) 
else 

lone = int(J0.0d0 + 4.0d0*zabs(arg)/3.0d0) 
qntone — l.OdO 
do 10 k = l,lone 

fmone = dble(lone - k + 1) 

qntone = l.OdO - qntone*(0.5d0*arg)**2/(fmone*(fmone+fm)) 
10 continue 

besjl = f*qntone 
endif 

return 

end 
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function cdfact(n,arg) 
c calculate (arg/2) A n/n! 
integer n,k 

double complex arg,ic, cdfact, f 
double precision theta, ertolp, ertolm 

common /DERRTOL/ ertolp,ertolm 

ic = (0.d0,l.d0) 
f = l.OdO 
if (n ,eq. 0) then 
cdfact = l.OdO 
else 

if (zabs(arg) .It. 1.5d0*ertolm) then 
cdfact = ertolm 
return 
endif 

do 5 k = l,n 

f = f*arg/(2.0d0*dble(n-k+l)) 
theta = datan(dimag(f)/dble(f)) 
if (zabs(f) .It. ertolm) then 

cdfact = ertolm*exp(ic*dbIe(n)*theta) 
return 
endif 

if (zabs(f) .gt. ertolp) then 

cdfact = ertolp*exp(ic*dble(n)*theta) 
return 
endif 

5 continue 
cdfact = f 
endif 

return 

end 

£* **** ****** * ** * * + * + * + + + * ** i:#****** * ******* * 

subroutine dcdbsjl(arg,m,dbesjl) 
integer m,lexit,k 

double complex arg,dbesj l,dqn ton, f, cdfact, Ic,besjl 
double precision Fin, amexit, theta, ertolp, ertolin 

common /DERRTOL/ ertolp, ertolm 


ic = (0.d0,l,d0) 
if (m .eq. 0) then 
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call dcbsjl(arg,l,besjl) 
dbesj 1 — -besj 1 
return 
endif 

f = cdract(m,arg) 
if (zabs(f) .It. 1.5dO*ertolm) then 
dbesj 1 = O.OdO 
return 
endif 

if (zabs(f) .gt. .5d0*ertolp) then 
theta = datan(dirnag(f)/dble(f)) 
dbesj 1 = ertolp*exp(ic*dble(m)* theta) 
return 
endif 

frn = dble(m) 

lexit = in t( IG.OdO + 1.6d0*zabs(arg)) 
dqnton = l.OdO 
do 10 k = 1, lexit 

amexit = dble(lexit-k+l) 

dqnton = l.OdO - dqnton*(0.5d0*arg)**2*(fm+2.0d0*amexit) / 
(amexit* (fm4*amexit)*(fm+2.0d0*amexit-2.0d0)) 

10 continue 

dbesj 1 = 0.5d0*dqnton*cdfact(rn-l,arg) 

return 

end 

c ******************************************* 
subroutine dcbsy l(arg,rn,besyl) 

integer m,k,lz,klz,]flag 

double complex arg,besyl,emz,besjl,pemzj,qntbz 
double complex cdfact,f,fmz,qntdk,sumdk,ic,gmz 

double precision fm,gamma,pi,repi,am,ertolp,ertolm,fklz,dk,ddigam 
double precision theta 

common /DERRTOL/ ertolp,ertolm 

ic = (0.d0,l.d0) 
f “ cdfact(m,arg) 
frn = dble(m) 

gamma = 0.57721566490153d0 
pi = 3.14l5926535897932d0 
repi = l.OdO/pi 

if (zabs(arg) .It. 1.5d0*ertolr») then 
besyl = -ertolp 
return 
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endif 


emz = 2.0d0*repi*(gamma + zlog(0.5d0*arg)) 
call dcbsjl(arg,m,besjl) 
pemzj = emz*besjl 

if (m .eq. 0) then 
fmz = O.OdO 
else if (ni .eq. 1) then 
fmz = -2.d0*repi/arg 
else if (zabs(f) .It. 1.5d0*ertolm) then 
theta — datan(dimag(f)/dble(f)) 
fmz = -erfcolp*exp(ic*fm*theta) 
else if (zabs(f) .gt. .5d0*ertolp) then 
fmz = O.OdO 
else 

qntbz = l.OdO 
do 10 k = l,m-i 
am = dble(m-k) 

qntbz = l.OdO + qntbz*(0.5d0*arg)**2/(am*(fm-ain)) 

10 continue 

fmz = -qntbz*repi / (f*fm) r i 
endif 

Iz — int(7.0d0+ 1.5d0*zabs(arg)) 
qntdk = l.OdO 
do 20 k = l,lz 
klz = lz - k -f 1 
fklz — dble(klz) 

dk = (ddigam(klz + 1) + ddigam(m+klz+l)) / 

((ddigam(klz)+ddigam(klz+m))*(fklz+l.dO)*(fmd-fklz+i.dO)) 
qntdk = l.OdO - qntdk*dk*(0.5d0*arg)**2 
20 continue 

suindk — ddigain(m)-(qntdk*(l.d0-fddigam(m+l))*(arg*0.5d0)**2 / 
(fm+ l.OdO)) 

if (zabs(f) .It. 1.5d0*ertolm) then 
gmz = O.OdO 

else if (zabs(f) .gt. .5d0*ertolp) then 
theta = datan(dimag(f)/dble(f)) 
gmz = ertolp*exp(ic*fm*theta) 
else 

gmz = -sumdk*repi*f 
endif 

besyl = pemzj + fmz + gmz 

call cdtolch (besy 1 ,besy 1 ,erfcolp, O.OdO, iflag) 

return 

end 
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c ***************** ************************** 
subroutine dcdbsy 1 (arg,m,dbesy 1) 

integer m,k,kmm,larg,klarg 

double complex arg,dbesy l,emz,dbesjl,demzjl,besjl 
double complex demzj2, qntdb, f, dfmz, qntder,sumder 
double complex dergmz,cdfact,ic 

double precision fm, gamma, pi, repi,fk,flarg,dkder,ddigam, ertolp 
double precision theta, ertolm 

common /DERRTOL/ ertolp, ertolm 

ic = (0.d0,l.d0) 
f = cdfact(m,arg) 
fm = dble(m) 

gamma = 0.57721566490153d0 
pi = 3.1415026535B97932d0 
repi = l.OdO/pi 

if (zabs(arg) .It. l.5*ertolrn) then 
dbesyl = ertolp 
return 
endif 

ernz = 2.0d0*repi* (gamma + zlog(0.5d0*arg)) 

call dcdbsjl(arg,m,dbesjl) 

dernzji = emz*dbesjl 

call dcbsj 1 (arg,m,besj 1 ) 

demzj2 = besji*2.0d0*repi/arg 

if (m .eq. 0) then 
dfmz = O.OdO 
else if (in .eq. 1) then 

dfmz = 2.0d0*repi/arg**2 
else 

if (abs(arg*f) .It. .5*ertolrh) then 
qntdb = ertolp 
else 

qntdb = 1 .d0/ (.5d0*arg*f) 
endif 
k = 1 

10 kxnm = m - 1 -.k 
fk = dble(k) ' iU ' 

qntdb = qntdb f ((fm-(2.d0*fk))*cdfact(k,arg)) / 
(cdfact(kmm,arg)*(.5d0*arg)**2) 
if (zabs(qntdb) .It. 1.5d0*ertolm) then 
dfmz = O.OdO 

else if (zabs(qntdb) .gt. .5d0*ertolp) then 
theta = datan(dimag(qntdb)/dble(qntdb)) 
dfmz = ertolp*exp(ic*fni*theta) 
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else 

dfmz = .5dO*qntdb*repi 
endif 
k = k + 1 

if (k .le. m-1) go to 10 
endif 

if (zabs(f) .It. 1.5d0*ertolm) then 
dergmz = O.OdO 

else if (zabs(f) .gt. ,5d0*ertolp) then 
theta = datan(dimag(f)/dble(f)) 
dergmz = ertolp*exp(ic*fm*theta) 
else 

larg = int(8.0d0+1.4d0*zabs(arg)) 
qntder = l.OdO 
do 20 k ~ l,Iarg 

klarg — larg - k + 1 
flarg = dble(larg-k + l) 

dkder = (ddigam(klarg+1) + ddigam(m-fklarg+ 1)) * 
(fm+2.d0*flarg+2.d0) 

. / ( (ddigam(klarg) fddigam(m-bklarg))* (flarg +I.d0)* 

(fm+flarg 4d.dO)*(fm-f2.dU*flarg) ) 
qntder = l.OdO - qntder*dkder*(0.5d0*arg)**2 
20 continue 

sumder = fm*ddigam(m) - qntder*(l.dO+ddigam(m+l))* 

(fm h 2.d0)*(0.5d0*arg)**2/(fm+l.d0) 
dergmz = -0.5d0*repi*f*sumder/(0.5d0*arg) 
endif 

dbesyl — demzjl + deinzj2 + dfmz + dergmz 
call cdtolch(dbesyl,dbesyl,ertolp,O.OdO,iflag) 
return 
end 

subroutine dcbesjy3(arg,m,icode,besj3,besy3,dbesj3,dbesy3) 
c Airy function approximation to Bessel function, complex argument 
integer m,icode,lrnk,k 

double complex arg,besj3,besy3,dbesj3,dbesy3,srarg,angle 
double complex qnta,qntb,uink,vmk,qntad,qntbd,rmk,smk 
double precision ertolp,ertoim,fm,pi,fmsr,flmk4,amk,bmk,arnkd,bmkd 
double precision theta 

common /DERRTOL/ ertolp,ertolm 

fm = dble(m) 

pi = 3.141 5926535897932d0 
srarg — sqrt(2.0d0/(pi*arg)) 
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angle = arg - 0.25d0*pi*(1.0d0 + 2.0d0*fm) 
fmsr = 4.0d0*fm**2 

link = int( 9.d0 + zabs(arg)*dexp(-.0458d0*zabs(arg)) ) 
if (Irak .It. 1) lmk = 1 

qnta = l.OdO 
qntb = l.OdO 
do 10 k = l,lmk 

{lmk4 — 4.0d0*dble(Jmk-k + 1) 

amk = (4.dO*(frnsr-(flmk4-3.dO)**2)*(fnisr-(flmk4-l.dO)**2)) / 
(flmk4*(flriik4 - 2.dU)) 

bmk - (4.d0*(fmsr-(flrnk4-l.d0)**2)*(fmsr-(flmk4-hl.d0)**2)) / 

. (flrnk4*(flmk4 + 2.d0)) 

qnta = l.dO - qnta*arnk/(8.d0*arg)**2 
qntb = l.dO - qntb*brnk/(8.d0*arg)**2 
if (zabs(qnta) .ge. .5d0*ertolp) then 

theta ~ datan(dimag(qnta)/dble(qnta)) 
qnta = ertolp + exp((0.d0,l.d0)*theta) 
endif 

if (zabs(qntb) .ge. .5d0*ertolp) then 

theta — datan(dimag(qntb)/dble(qntb)) 
qntb — ertolp*exp((0.d0,l.d0)*theta) 
endif 

10 continue 
umk — qnta 

vmk — qntb*(fmsr - 1.0d0)/(8.0d0*arg) 

beaj3 = sra.rg*(urnk*cos(angle) - vmk*sin(angle)) 

if (zabs(besj3) .ge. .5d0 4: ertolp) then 

theta = datan(dimag(besj3)/dble(besj3)) 
besj3 = ertolp*exp((0.d0,l.d0)*theta) 
endif 

besy3= srarg*(urnk*sin(angle) + vmk*cos(angle)) 
if (zabs(besy3) .ge. ,5d0*ertolp) then 

theta = datan(dimag(besy3)/dble(besy3)) 
besy3 = ertolp*exp((0.d0,l.d0)*theta) 
endif 

if (icode .eq. 3) (lien 
qntad = l.OdO 
qntbd = l.OdO 
do 20 k = ljlmk 

flmk4 = 4.0d0*dble(lmk-k+l) 

amkd = (4.d0*(finsr+(flrnk4**2-l.d0))*(fmsr“(flmk4-3.d0)**2)* 
(finsr-(flmk4-5.d0)**2)) / 

(flmk4*(flmk4 - 2.d0)*(fmsr+(flmk4-4.d0)**2-l.d0)) 
brnkd = (4.d0*(fmsr+(flmk4-j-2.d0)**2-l.d0) + 

(fmsr-(flmk4-3.d0)**2)*(friisr-(flmk4-l.d0)**2)) / 
(flrnk4*(flmk4 + 2.d0)*(fmsr+(flmk4-2.d0)**2 - l.dO)) 
qntad = l.dO - qntad*amkd/(8.d0*arg)**2 
qntbd = l.dO - qntbd*bmkd/(8.d0*arg)**2 
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if (zabs(qntad) .ge. .5dO*ertolp) then 

theta = datan(dimag(qnlad)/dble(qntad)) 
qntad = erfcolp*exp((0.d0,l.d0)*theta) 
endif 

if (zabs(qntbd) .ge. .5d0*ertolp) then 

theta = datan(dirnag(qntbd)/dble(qntbd)) 
qntbd = ertoIp*exp((0.dO,l.dU)*theta) 
endif 

20 continue 

rink — qntad 

snik = qntbd*(fmsr + 3.0d0)/(8.0d0*arg) 
dbesj3 = -srarg*(rmk*sin(angle) + smk*cos(angle)) 
if (zabs(dbesj3) .ge. .5dU*ertolp) then 

theta = datan(dimag(dbesj3)/dble(dbesj3)) 
dbesj3 = ertolp*exp((0.d0,l.d0)*theta) 
endif 

dbesy3 = srarg*(rmk*cos(angle) - smk*sin(angle)) 
if (zabs(dbesy3) .ge. .5d0*ertolp) then 

theta = datan(dimag(dbesy3)/dble(dbesy3)) 
dbesy3 — erto!p*exp((0.d0,l.d0)*theta) 
endif 
endif 

return 

end 

subroutine dcbesjy2(arg,m,icode,besj2,besy2,dbesj2,dbesy2) 

c asymptotic approximation to Bessel function, complex argument 
c Abr-Stegun 9.3.35, 9.3.36 

integer m,icode 

double complex arg,besj2,besy2,dbesj2,dbesy2,x,zeta,ql 
double complex qli,g3,g32,g2,gsq,al,bO,cO,dl,h,y,airy 
double complex d,ccos, phi, fi,f2,f3,biry, dairy ,dbiry,ic 
double complex ul,u2,u3,bl,qq,u4,u5,a2,b2 
double precision ertolp,ertolm,fm,q,w,a,b 

common /DERRTOL/ ertolp,ertolm 

ic = (0.d0,l.d0) 
fm — dble(m) 
x = arg/fm 

if (zabs(x-l.dO) .It. 1.5d0*ertolm) then 
zeta — O.OdO 

else if (dble(x) .gt. l.dO) then 

q — dble(x) / (dble(x)**2H-dimag(x)**2) 
w = -dimag(x) / (dble(x)**2+dimag(x)**2) 
c q + iw = 1/x 
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a = ,5dO*( sqrt((q+l.dO)**2+w**2) + sqrt((q-l.dO)**2+w**2) ) 
b = .5dO*( sqrl((q } l.dO)**2-fw**2) - sqrt((q-l.d0)**2+w**2) ) 
d = a + sqrt(a**2-l.dO) 
ccos = acos(b) - ic*zlog(d) 

if (dimag(x) .gt. O.dO) ccos = dble(ccos)-ic*dimag(ccos) 
zeta = -(1.5d0*(sqrt(x**2-l.d0) - ccos))**(2.dO/3.dO) 
else 

zeta = (1.5d0*(zlog((l .OdO +sqrt(l.d0-x**2))/x) 

- sqrt(l.dO-x**2)))**(2.dO/3.dO) 

endif 

ir (zabs(x-l.dO) .It. 0.02d0) then 
h = l.OdO - x 

phi = 2.0d0**(l.d0/3.d0)*(1.0d0 + 0.2d0*h + 3.d0*h**2/35.d0 
+ 73.d0*h**3/1575.d0 + 35209.d0*h**4/1212750.d0 + 
380069.d0*h**5/18768750.d0) 
al = -1. dO/225. dO - 71.d0*h/38500.d0 
bO = (1. d0/70. dO + 2.d0*h/225.d0)*2**(l.d0/3.d0) 
bl = 2.d0**( t.d0/3.)*(-l 2l3.d0/1023750.d0 - 3757.d0*h/2695000. 

- h**2*(8.9962899979797d-4 + h*(.0002753433716d0 - h* 
(,00018048868d0 -|- h*.0004108523)))) 
a2 = 6.937355413546877d-4 + h* (.00046448349036601 - li* 
(,0002890362546053d0 + h*(.0008747649439535d0 + 
h*.00102971637614))) 

b2 = -2.d0**(l./3.)*(4.3829l8094497229d-4 + h* 

(7. 11048651 1691 ld-4 + h*5.318984348085d-4)) 
b2 = b2 - 2.d0**(l./3.)*h**3*2.t82958472d-4 
cO = (O.ldO + 0.02d0*h)*2**(2.d0/3.d0) 
dl = 23.d0/3150,d0 + 1453.d0*h/346500.d0 
else 

ql = zeta/ ( 1 .0d0-x**J2) 
phi = (4.0d0*ql)**0.25d0 • 

fl = l.OdO - x**2 
f2 = fl**2 
f3 = 11**3 
qli = l.OdO/ql 
g3 = qli**3 
g32 = qli**1.5d0 
g2 = qli**2 
gsq = sqrt (qli) 

al = (-455.d0*g3/4608.d0-7.d0*g32*(fl-5.d0/3.d0)/384.d0+385.d0/ 
1152. dO - 77.dO*fl/192.dO + 9.d0*f2/128.d0 ) / f3 
bO = (-5.dO*g2/48.dO + gsq*(5.dO/24.dO - fl/8.dO))/f2 
cO = (7.dO*qli/48.dO - 3.dO*sqrt(ql)*(7.dO/72.dO - fl/8.dO))/fl 
dl = (385.d0*g3/4608.d0+5.d0*g32*(-fl+7.d0/9.d0)/128.d0 
-15.dO*f2/128.dO + 33.dO*fl/64.dO - 455.d0/1152.d0 ) / f3 
qq = sqrt(fl/zeta) 

ul = qq*(i.dO-5.dO/(fl*3.dO)) / (8.dO*fl) 
u2 = (4.5dO - 77.dO/(fl*3.dO)*(l.dO-5.dO/(fl*6.dO))) / 

(64.d0*fl) 


57 



u3 = qq*(75.d0/2. - 456.3d0/fl + 17017.d0/(f2*18.d0)*(l.d0- 
5.dO/(fl*9.dO))) / (512.dO*P2) 

u4 = (3675.dO/8.dO - 9683.3d0/fl + 2717.dO/f2*(53.dO/4.dO - 
2737.d0/(fl*162.d0)*(l.d0-5.d0/(12.d0*fl))))/(4096.d0*f2) 
u5 = 59535. dO/8. dO - 221.d0/(4.d0*fl)*(305923.d0/70.d0- 
77.dO/fl*(14743.dO/45.dO - 95.d0/fl*(67.d0/9.d0- 
3335. dO/ (486.dO*fl)*(l.dO-l.dO/ (3.dO*fl))))) 
u5 = qq*u5/(f3*32768.dO) 

bl = -u3 - 5.d0/(zeta**2*48.d0)*(u2+77.d0/(zeta*96.d0)* 

(ul + 221.dO/(zeta**2*144.dO))) 
b2 =-u5 - 5. dO/ (zeta**2*48.dO)*(u4 4 - 77.d0/(zeta*96.d0)* 

(u3 + 221. dO/ (zeta**2*144.dO)*(u2 + 437.dO/(zeta*192.dO)* 
(ul -I- 145. dO/ (zeta**2*48,dO))))) 
a2 = u4 - 7.dO/(zeta*48.dO)*(u3 + 65.d0/(zeta**2*96.d0)*(u2 + 
209.d0/(zeta*144.d0)*(ul + 425.dO/(zeta**2*192.dO)))) 

endif 

y = zeta*fm**(2.dO/3.dO) 

call cdairy fn(y, airy ,biry, dairy ,dbiry) 

if (icode .eq. 1) then 

besj2 = phi/(fm**(l.dO/3.dO))*(airy*(l.dO-fal/fm**2+a2/fm**4) 
+ dairy / (fm**(4.dO/3.dO)J*(bO+bl/fm**2-fb2/fm**4)) 
else if (icode .eq. 2) then 

besy2 — -phi/(fm**(l.dO/3.dO))*(biry*(l.dO+al/frn**2+a2/fm**4) 
+ dbiry / (fm**(4.dO/3.dO))*(bO+bl/fm**2+b2/fm**4)) 
call cdtolch(besy2,besy2,ertolp,0.0d0,iflag) 
else 

besj2 = phi/(fm**(l.d0/3.d0))*(airy*(l.d0 fal/fm**2-f a2/fm**4) 
+ dairy / (fni**(4.dO/3.dO))*(bO-fbl/fm**2-l-b2/fm**4)) 
besy2 - - P hi/(fm**(l.d0/3.d0))*(biry*(l.d0+al/fm**2+a2/fm**4) 
+ dbiry / (fm**(4.dO/3.dO))*(bO+bl/fm**2+b2/fm**4)) 
call cdtolch(besy2,besy2,ertolp,0.0d0,iflag) 
dbesj2 =-2.d0/( fm**(2.dO/3.dO)*x*phi)*(airy/ 

(frri**(2.d0/3.d0))*c0 + dairy*(l.dO + dl/(fm**2.d0))) 
dbesy2 = 2.d0/(fm**(2.d0/3.d0)*x*phi)*(biry*c0/ 

(fm**(2.dO/3.dO)) + dbiry*(l.dO + dl/(fm**2.d0))) 
call cdtolch(dbesy2,dbesy2,ertolp,0.0d0,iflag) 
endif 


return 

end 

subroutine cdairyfn(y,airy,biry, dairy, dbiry) 

integer icode,ll,k,k2,l3,l4,l5,lp,lpt 
double complex y, airy, biry, dairy, dbiry, flam,reflam 
double com|)lex clam,srtwo,qntv,ffmi,vffml,vffmr,qntw 
double complex qntdv,ffin3,dvffmr,qntdw,ffm4,dwffmr,suml 
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double complex sum4,diftl,sumtl,dift2,sumt2, qntf, qntg 
double complex ff,fg,fdf,fdg,fml,fpl,finp,gml,gpl,fmpt 
double complex al,bl,dal,dbl,a2,b2,da2,db2,a3,b3,da3,db3 
double complex aby, slam, ffm2,wffmr,sum2,sum3, qntdf 
double complex qntdg, sxfmpt,di2,dl3,sxfmp 
double precision ertolp,ertolm,cl,c2,ffm5,err,rl2,rl3,pi,srpi 

common /DERRTOL/ ertolp, ertolm 

cl = 0.355028053887817d0 
c2 = 0.258819403792807d0 

if (zabs(y) .It. ertolm) then 

c Ai(0) etc from Abr-Stegun 10.4.4, 10.4.5 
airy — cl 

biry = cl*sqrt(3.d0) 

dairy = -c2 

dbiry = sqrt(3.d0)*c2 

else if (zabs(y) .ge. 25. do) then 

c at machine tolerance 

airy — • O.OdO 
dairy = O.OdO 
hiry = ertolp 
dbiry = ertolp 

else if (zabs(y) .It. 4.0d0) then 

c power series 

L5 = int(5.0d0 + 3.5d0*zabs(y)) 
qntf = l.OdO 
qntg = l.OdO 
qntdf — l.OdO 
qntdg — l.OdO 
do 84 k = 1,L5 
ffm5 = 3.dO*dble(L5-k+l) 
qntf = l.OdO + qntf*y**3/(ffm5*(ffm5 - l.OdO)) 
qntg = l.OdO + qntg*y**3/(ffm5*(ffm5 + l.OdO)) 
qntdf = l.OdO + qnfcdf*y**3/(ffm5*(ffm5 + 2.0d0)) 
qntdg = l.OdO + qntdg*y**3/(ffm5*(ffm5-2.0d0)) 

84 continue 
ff = qntf 
fg = qntg*y 
fdf = qntdf*0.5dl)*y**2 
fdg = qntdg 
airy = cl*ff - c2*fg 
dairy = cl*fdf - c2*fdg 
biry = sqrt(3.0d0)*(cl*ff + c2*fg) 
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dbiry = sqrt(3.0d0)*(cl*fdf + c2*fdg) 
call cdtolch(airy,airy,ertolp,O.OdO,iflag) 
call cdtolch(biry,biry,ertolp,O.OdO,iflag) 
call cdtolch (dairy, dairy ,ertolp,O.Od0,iflag) 
call cdtolch(dbiry, dbiry ,ertolp,0.0d0,iflag) 

else if ( (dble(y) .ge. O.OdO) .or. 

(dabs(dimag(y)/dble(y)) .ge. sqrt(3.dO)/2.dO) ) then 

exponential approx (Abr-Stegun 10.4.59) 

flam = (2.d0/3.d0)*(y)**1.5d0 

reflam = l.OdO/flam 

pi = 3.141592G535897932d0 

srpi = dsqrl(pi) 

aby = (y)**0.25d0 

lp = int(5.d0 + 70.0d0*zabs(reflam)) 

fml = l.OdO 

fpl = l.OdO 

do 115 k = l,lp 

sxfrnp = 6.0d0*dble(lp-k-f-l) 

fml — I.d0-fml*reflam*(sxfinp-5.d0)*(sxfmp-l.d0)/ 
(12,d0*sxfmp) 

fpl = I.d0+fpl*reflam*(sxfmp-5.d0)*(sxfmp-l.d0)/ 

. (12.d0*sxfmp) 

continue 

Ipt = int(4,0d0 + 90.0d0*zabs(reflam)) 
gml — l.OdO 
gpl = l.OdO 
do 125 k = l,lpt 

sxfmpt = 6.0d0*dble(lpt-k + l) 

gml = I.d0-gml t reflam + (sxfmpt+l.d0) i£ (sxfmpt-7.d0)/ 
(22.d0*sxfrnpt) 

gpl ~ LdO+gpl*refWm*(sxfmpt-bl.dO)*(sxfmpt-7.dO)/ 
(12.d0*sxfmpt) 
continue 

airy = 0.5d0*y**(-0.25d0)*fml*exp(-flam)/srpi 
dairy = -0.5d0*y**0.25d0*gml*exp(-flam)/srpi 
biry = y* :f: (-0.25d0)*fpl*exp(flam)/srpi 
dbiry =y**0.25d0*gpl*exp(flam)/srpi 
call cdtolch (airy, airy ,ertolp,0.0d0,iflag) 
call cdtolch(biry,biry,ertolp,0.0d0,iflag) 
call cdtolch (dairy, dairy ,ertolp,0.0d0,iflag) 
call cdtolch(dbiry,dbiry,ertolp,O.OdO,iflag) 

else 

trig approx (Abr-Stegun 10.4.60) 
flam = (2.dO/3.dO)*(-y)**1.5dO 
reflam — l.OdO/flam 
aby = (-y)**0.25d0 



slam = sin(flarn) 

clam = cos(flam) 

pi = 3.141 5926535897932d0 

srpi = dsqrt(pi) 

srtwo = sqrt(2.dO)/2.dO 

LI = int(1.0d0 + 75.0d0*zabs(reflam)) 

qntv = l.OdO 
do 20 k = 1,L1 

ffinl = 12. dO *dble(L l-k+1) 

vffmr = (ffinl- 1 1 ,dO)*(ffml-7.dO)*(ffml-5.dO)*(ffml-l.dO) 

/ (144.d0*ffml*(ffml-6.d0)) 
qntv = l.OdO - (qntv*vffmr*reflam**2.d0) 

20 continue 

L2 = int(J.0d0 + 70.0d0*zabs(reflam)) 
qntw = l.OdO 
do 30 k = 1,L2 

ffm2 = 12.dO*dble(L2-k+l) 

wffmr = (ffrn2-5.d0) * ( Ffm2- 1 .dO) * (ffm2+ 1 .dO)* (ffm2-|-5.d0) 

/ ( 1 44.d0*ffm2*(ffm2 + G.dO)) 
qntw = l.OdO - (qntw*wffmr*reflam**2.d0) 

30 continue 

L3 = int(1.0d0 + 60.0d0*zabs(reflam)) 
qntdv — l.OdO 
do 40 k = 1,L3 
ffm3 = 12.dO*dble(L3-k+l) 

dvffinr = (ffm3-7.d0) ,: (ffm3-l.dO)*(ffm3+l.dO)*(ffrn3+7.dO) 
/ (144.d0*ffm3*(ffm3 + 6.d0)) 
qntdv = l.OdO - (qntdv*dvffmr*reflam**2.d0) 

40 continue 

L4 = inl(1.0d0 + 33.0d0*zabs(reflam)) 
qntdw = l.OdO 
do 50 k = 1,L4 

ffm4 = 12.dO*dble(L4-k + l) 

dwffmr = (ffm4-13.dO)*(ffm4-7.dO)*(ffm4-5.dO)*(ffm4+l.dO) 
/ (144.d0*ffm4*(ffm4 - 6.d0)) 
qntdw = l.OdO - (qntdw*dwffmr*reflam**2.d0) 

50 continue 

suml = qntv /srpi 

sum2 — 5.d0*qntw*reflam/(72.d0*srpi) 

sum3 = -7.dO*qntdv*reflam/(72.dO*srpi) 

suin4 = qntdw/srpi 

diftl — srtwo*(suml - sum2) 

sumtl = srtwo*(suml + sum2) 

dift2 = srtwo*(sum4 - suin3) 

sumt2 = srtwo* (sum4fsuin3) 

airy — (diftl*clarn + surntl*slam)/aby 
dairy = (dift2*slam - sumt2*clam)*aby 
biry = (sumtl*clam - diftl*slam)/aby 
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dbiry = (sumt2*slam + dift2*clam)*aby 
call cdtolch(airy , airy ,ertolp, O.OdO, iflag) 
call cdtolch(biry,biry,ertolp, O.OdO, iflag) 
call cdtolch (dairy, dairy ,erto!p, O.OdO, iflag) 
call cdtolch(dbiry, dbiry ,ertoIp, O.OdO, iflag) 

endif 

return 

end 

subroutine cdtolch(cin,cout,setbig,setsmaIl, iflag) 
c complex, double precision check to avoid over- and underflow 
integer iflag 

double precision setbig,setsmall, theta, pi, ertolp,ertolm 
double complex cin,cout,ic 

common /DERRTOL/ ertolp,ertolm 

pi = 3.1415926535897932d0 
ic = (0.d0,l.d0) 
iflag = 0 

if (zabs(cin) .It. 1.5d0*ertolm) then 
theta = datan(dirnag(cin)/dble(cin)) 
if (dble(cin) .It, O.OdO) theta = theta + pi 
coat = setsmall*exp(ic*theta) 
iflag = 1 

else if (zabs(cin) ,gt. ,5d0*ertolp) then 
theta = datan(dimag(cin)/dble(cin)) 
if (dble(cin) .It. O.OdO) theta = theta + pi 
cout = setbig*exp(ic*theta) 
iflag = 1 
endif 
return 
end 
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TABLE I - SINGLE PRECISION, REAL ARGUMENT 


m 

X 


Y„(x) 

J-M 

viW 

1 

10.0 

4.3472745xl0* 2 

0.2490154 

-0.2502830 

3.0769626xl0' 2 

10 

2.0 

2.5153864xl0* 7 

-129184.5 

1.2346503x10'* 

631362.9 

20 

25.0 

5.1994048xl0' 2 

0.1980408 

-0.1230286 

2.1158613xl0’ 2 


TABLE n - DOUBLE PRECISION, REAL ARGUMENT 


m 

X 


Y m (x) 



Ym(x) 

1 

10.0 

4.3472746168805587xl0' 2 

0.2490154242067848 

-0.2502830390682086 

3.0769624863030031xl0‘ 2 

10 

2.0 

2.5153862827167365xl0’ 7 

-129184.5422080393 

1.2346503x10'* 

1.2346502937746957x10'* 

20 

25.0 

5. 19940492 283 02969x 1 0‘ 2 

0.1980407477628901 

-0.123-285643023131 

2.11586141 185 14424xl0' 2 


TABLE IB - SINGLE PRECISION, COMPLEX ARGUMENT 


m 

r 

degrees 

JJx) 

Y m (x) 

CM 

Y» 

1 

10.0 

37 

48.46594, 14.90409 

-14.90460,48.46558 

12.67906, -47.47811 

47.47851.12.67855 

10 

2.0 

-15 

-2.1449186xl0' 7 , -1.3728736xl0' 7 

1064829.7, -69704.96 

-8.3601503xl0' 7 , -9.3365060x1 O' 7 

-594765.0, 188178.8 

n 

20 

25.0 

80 

1309767., 1275952. 

-4.3071617xl0' 9 , 3.3826593xl0' 9 

1714659., -1538981. 

-3. 9587 404x1 0' 9 , 5,8104996x10 9 


TABLE IV - DOUBLE PRECISION, COMPLEX ARGUMENT 


m 

r 

*> 

degrees 

J m (x) 

Y m (x) 

1 

10.0 

37 

48.46594122869320, 14.90408497152226 

-14.90459530350984, 48.46557450542743 

10 

2.0 

-15 

-2.1449186196047720xl0’ 7 , -1.3728736108270844xl0' 7 

-8. 360 15026034 10100x1 O' 7 , -9.336506171 1294054xl0' 7 

20 

25.0 

80 

1309766.766321728,1275951.977909705 

1714659.236311793, -1538981.266957863 


Double Precision, Complex Argument 


m 

r 

degrees 

Cm 

Y'(x) 

B 

10.0 

37 

12.67906150921442, -47.47811148316850 

47.47851111589829, 12.67855209985360 


2.0 

-15 

-8.30615026034101000xl0* 7 , -9.3365O6l711294057xl0' 7 

-594764.9878571380,188178.8391119776 


25.0 

80 

1714659.236311793, -1538981.266957863 

-3.958740335483453 lxlO' 9 , 5.8104997303197833xl0* 9 
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Figure 1 —Evaluation regions for Bessel functions with real 
argument. 


Figure 2— Evaluation regions for 
Bessel functions with complex 
argument. 


Trigonometric 



Exponential 

expansion 

t 

Power series 

I 

expansion 

1 

-10 

I 

0 

1 

6 



Figure 3.— Evaluation regions for Airy functions with real 
argument. 



Figure 4. — Evaluation regions for Airy functions with 
complex argument. 
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